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CHAPTER I 


Introduction 


Microstrip patch antennas offer considerable advantages in terms of weight, aero- 
dynamic drag, cost, flexibility and observables over more conventional protruding 
antennas. These flat patch antennas were first proposed over thirty years ago by 
Deschamps [1] in the United States and Gutton and Baisinot [2] in France. Such 
antennas have been analyzed and developed for planar as well as curved platforms. 
However, the methods used in these designs employ gross approximations, suffer 
from extreme computational burden or require expensive physical experiments. The 
goal of this thesis is to develop accurate and efficient numerical modeling techniques 
which represents actual antenna structures mounted on curved surfaces with a high 
degree of fidelity. 

1.1 Background 

Microstrip patch antennas are conformal metallic elements printed on a grounded 
dielectric substrate. Although a single patch element offers low gain as compared to 
other antennas, an array of such elements offers considerable gain and pattern shaping 
flexibility in a low profile package. Traditionally, such antennas were designed using 
either expensive measurements or approximate formulae such as the cavity model. 
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The cavity model [3] represents the radiating structure by magnetic current walls 
placed around the perimeter of the patch between the groundplane and the surface 
of the substrate. The original cavity model assumed a planar-rect angular patch but 
it has been extended to cylindrical-rectangular patches [4. 5]. The cavity model has 
proven quite accurate for radiation pattern calculations near resonance. However, 
it is not suitable for operation away from resonance since it is based upon a single 
mode assumption. Additionally, it is not especially accurate for input impedance 
calculations since probe feeds typically excite multiple modes within the cavity. The 
cavity model also ignores mutual and surface wave coupling between antenna ele- 
ments which is often a concern to array designers since it results in parasitic energy 
loss. Finally, the cavity model cannot be used for scattering calculations. 

To overcome the deficiencies of the cavity model, rigorous integral equation for- 
mulations have been investigated [6]- [9]. Integral equations utilize an appropriate 
Green’s function and equivalent currents to represent the fields. Such an approach 
offers a rigorous representation of the local fields which includes mutual coupling and 
surface wave effects. However, this approach suffers from two primary drawbacks: 
computational expense and memory demand. In particular, accurate computation of 
the substrate Green’s function is both difficult and expensive since an accurate eval- 
uation of a Sommerfeld integral is required [10]. Nevertheless, the integral equation 
technique has been extended to patch antennas printed on a coated circular cylin- 
der [11, 12], which requires an extremely complex substrate dyadic Green’s function. 
To reduce the memory demand and matrix fill time from 0(N 2 ) to O(N), the Bi- 
conjugate Gradient-Fast Fourier Transform (BiCG-FFT) may be used provided the 
substrate is either flat or circular and uniform discretization is used. 

Although integral equation formulations have proven quite useful for design and 
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analysis of patch antennas, they require an assumption which limits their fidelity for 
real-world antennas. Since a substrate Green s function is used in the formulation, 
the substrate and underlying metallic surface must be infinite. In practice, antenna 
designers place each radiator within a metallic cavity for suppression of unwanted 
coupling [13]. These structures have finite substrates which require a substantially 
more complex (and expensive) dyadic Green's function for their implementation. 
Thus, integral equation formulations cannot efficiently model such arrays and hence 
are of limited value in their design. To overcome the computational burden of integral 
equations, an alternative methodology was first proposed separately by Silvester and 
Hsieh [14] and McDonald and Wexler [15] which is now known as the Finite Element- 
Boundary Integral (FE-BI) method. Jin and Volakis [16] adapted the FE-BI tech- 
nique for planar cavity-backed patch antennas. In this approach, planar-rectangular 
cavities were discretized using brick-shaped finite elements and the computational 
domain was closed with a boundary integral which provided an exact mesh closure 
condition. This method was found to be capable of modeling large complex arrays 
of cavity-backed patch antennas with low memory and computational burden. In 
addition to radiation by finite, possibly aperiodic, arrays of such structures, Jin and 
Volakis also investigated their scattering behavior [17]. 

The FE-BI method has been applied to two-dimensional [18], Body of Revolution 
[19] and planar three-dimensional problems [16, 20]. The primary reason for the lack 
of interest in applying the FE-BI approach to more general curved geometries is that 
the method is no longer efficient in terms of memory or computational demand since 
FFTs may not be used in conjunction with the boundary integral. However, this 
technique can be applied to cylindrical bodies as well as planar geometries without 
excessive computational cost. A formulation can be developed which yields compa- 
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rable performance by utilizing special finite elements, dyadic Green's functions and 
the BiCG-FFT solver. 

The FE-BI technique is both accurate and efficient for cavities embedded in a 
ground plane or in a metallic cylinder. However, for any other platform, the com- 
putational burden of the FE-BI approach is excessive since the BiCG-FFT solver 
ma\ not be used. In fact, the only dyadic Green's function available for arbitrarv 
surfaces is the free-space Green’s function which would require both electric and 
magnetic currents over the entire mesh boundary. Accordingly, Absorbing Bound- 
ary Conditions (ABCs) have been proposed to provide an approximate mesh closure 
condition which avoids the use of an expensive boundary integral. ABCs have been 
used extensively for large scattering applications [21] but not for antennas. Unfortu- 
nately, most ABCs assume a spherical mesh boundary [22, 23] which would result in 
excessive unknowns. Chatterjee and Volakis [24, 25] have proposed conformal vec- 
tor ABCs which follow the shape of the underlying body and in this work they are 
applied to scattering and radiation by cavity-backed conformal antennas which are 
recessed in a metallic cylinder. Conformal ABCs have also recently been proposed 
by Stupfel [26]. Verification of the accuracy of the FE-ABC method for antenna ap- 
plications is important since it is the only practical means of modeling large doubly 
curved and possibly coated conformal arrays. 

1.2 Overview 

In this thesis, the FE method is extended to cavity-backed conformal antenna 
arrays embedded in a circular, metallic, infinite cylinder. Both the boundary integral 
and absorbing boundary mesh closure conditions will be used for terminating the 
mesh. These two approaches will be contrasted and used to study the scattering and 
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radiation behavior of several useful antenna configurations. An important feature of 
this study will be to examine the effect of curvature and cavity size on the scattering 
and radiation properties of wraparound conformal antenna arrays. 

Chapter 2 introduces the fundamental electromagnetic principles used in the de- 
velopment of the FE-BI and FE-ABC formulations. Included are Maxwell's equa- 
tions, natural boundary conditions, the radiation condition, integral representation 
of fields in terms of the equivalence principal and the definition of radar cross sec- 
tion. In addition, the notation used throughout this thesis will be presented. An 
introduction to numerical formulations will be presented along wdth a discussion of 
matrix assembly. In particular, the integral equation and finite element methods will 
be compared with respect to matrix cost. 

In Chapter 3 the FE-BI formulation is developed for circular cylinders. Of par- 
ticular interest are the newly introduced cylindrical shell elements and an efficient 
evaluation of the cylindrical dyadic Green’s function. Appropriate plane wave, inte- 
rior source, far-zone field and input impedance formulae are given. The emphasis of 
this formulation is on efficiency in terms of both memory and computational demand. 
In view of this goal, the formulation will assume uniform surface gridding so that the 
BiCG-FFT solver may be used. Also, large radius asymptotic formulae for the cylin- 
drical dyadic Green’s function will be presented along with the traditional modal 
formulae and criteria will be established for the appropriate use of these formulae. 

Chapter 4 introduces the conformal ABCs for mesh termination. The FE-ABC 
equations are developed which utilize a novel mixed total/ scattered field formulation. 
This approach minimizes the potential error propagation usually associated with 
a total field formulation applied in the exterior region while retaining the lower 
computational burden of a total field formulation when used in the interior region. 
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New conformal vector ABCs which were developed by Ohatterjee and \olakis of 
the Radiation Laboratory are used for mesh closure. Both first and second order 
ABCs are presented for cylindrical-rectangular surfaces. The first order ABC always 
produces a symmetric system while the second order ABC is asymmetric. The second 
order ABC becomes symmetric for spherical boundaries and in the limit as the radius 
of the ABC surface becomes large. 

In Chapter 5 we investigate the scattering behavior of conformal antenna ar- 
rays. The FE-BI method is used to study the effect of curvature on the scatter- 
ing pattern and the resonant frequency. Both principal polarizations are discussed. 
Also, the effect of cavity size is examined. One array of particular interest is the 
wraparound array. The difference in scattering behavior between discrete and con- 
tinuous wraparound arrays is examined and related to known creeping wave and 
diffraction phenomena. The FE-ABC method will be compared to the FE-BI formu- 
lation and in particular, the minimum ABC displacement from the cavity aperture 
is determined. Such information is critical in minimizing the computational domain. 
Also, the FE-ABC method will be shown to recover the FE-BI data for monostatic 
scattering even when the transmitter/receiver is located behind the cylinder. 

Chapter 6 provides a corresponding study of radiation characteristics for con- 
formal patch antennas. The effect of curvature and cavity size will be determined 
for both axially and circumferentially polarized antenna elements. Additionally, the 
effect of curvature on the input impedance of both types of patches will be examined. 
An example of pattern coverage will be given for discrete, continuous and collar-type 
wraparound antennas. The pattern behavior as the number of antenna elements in- 
creases is demonstrated. Finally, the FE-ABC method is compared with the FE-BI 
method in order to determine the minimum required ABC displacement for accurate 
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radiation pattern calculations. The performance of the FE-ABC method for input 
impedance computations is also discussed. 



CHAPTER II 


Basic Electromagnetics and Numerical Techniques 

2.1 Basis Electromagnetics 

James Clerk Maxwell presented a set of six, coupled, time-dependent equations 
which completely describe all electromagnetic phenomena. Although Maxwell’s equa- 
tions allow for very general field variations in both space and time, in this work only 
time-harmonic fields will be considered where a e +Jwt time-dependency is assumed 
and suppressed. In addition, in this work only fields in a linear, isotropic medium 
will be considered. Accordingly, the following constitutive relations 

D = tE (2.1) 

B = fiH (2.2) 

have been used to develop Maxwell’s equations in differential or point form 

Vx£ = -Mi-jkZH (2.3) 

V x // = Ji+jkYE (2.4) 

where D is the electric flux density, E is the electric field, B is the magnetic flux 
density, H is the magnetic field, J, is an impressed electric current and M, is a 
fictitious magnetic current which is included so that (2.3) and (2.4) are symmet- 
ric. The medium is completely described by its intrinsic impedance (admittance) 
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Z = (yr) = yjf. where c and p denote the medium's permittivity and permeabil- 
ity, respectively. The wavenumber is denoted by k — = u where A is the 

wavelength and u : is the corresponding angular frequency. An important medium is 
free-space where k = k 0 = ^ with A 0 being the free-space wave length. 

Co = 8.854 x 10 ~ u F/m, p 0 = 47r x 10 H/m and the intrinsic impedance (admit- 
tance) is now denoted as Z 0 = (prj = ~ 377ft. 

Faraday’s (2.3) and Ampere-Maxwell’s (2.4) Laws are not sufficient for the com- 
plete specification of a general field. One of four dependent scalar equations must 
also be chosen w r hich may be derived from (2.3) and (2.4) 


v • E = p e /e 

(2.5) 


v • H = p m /p 

(2.6) 

V • J, = -jujpe 

(2.7) 


V ■ Mi = ju>p m 

. . - 

(2.8) 


In these, p e is the electric charge density (^) and likewise p m denotes the fictitious 
magnetic charge density (^r). The latter is not a physical quantity but its introduc- 
tion is important for mathematical constructions. Evidently, (2.5) and (2.6) are a 
statement of Gauss’ Law in differential form whereas (2.7) and (2.8) are referred to 
as the continuity or charge conservation equations. Typically, (2.3), (2.4) and either 

(2.5) or (2.7) are used to determine the field quantities. Such a solution is, however, 
not unique and will require the specification of additional constraints. 

Application of Maxwell’s equations (2.3) and (2.4) along with a derived equation 

(2.5) -(2.8) will result in an infinite number of solutions. A unique solution is obtained 
only after the imposition of various conditions such as boundary conditions, radia- 
tion conditions, edge conditions and transition conditions. The natural boundary 
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conditions satisfied by E and H at any dielectric interface are [27] 

hx{E x -E 2 ) = 0 (2.9) 

n x (//, - H 2 ) = 0 (2.10) 

where the subscripts denote the different media which are joined at the interface 
and n is the normal unit vector directed into medium one as shown in Figure 2.1. 
Medium one is the external medium while medium two is the internal medium. 

medium 1 


A 

n 



Figure 2.1: Material interface. 

The first condition (2.9) states that the tangential components of the electric field 
are continuous across the dielectric boundary while the tangential magnetic field 
components are also continuous for a dielectric. Boundary conditions may also be 
written in terms of normal field components 

n • — £ 2 ^ 2 ) = 0 (2.11) 

fi ■ (niHi - H2H2) = 


0 


( 2 . 12 ) 
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which simply states that the normal components of both the electric and magnetic 
fields are discontinuous if there is a discontinuity in the permittivity and permeabil- 
ity. respectively. Another important interface involves a dielectric (medium 1) and 
perfect electric conductor (medium 2). The tangential electric field must now vanish 
at such an interface 


h x E i=0 


(2.13) 


since a perfect electric conductor cannot sustain an internal field. Note that a surface 
current may now flow along the interface 

n x H] = J a (2-14) 


Similar conditions for the normal field components now may support a surface charge 
density. For perfect magnetic conductors, appropriate boundary conditions may be 
determined via duality. For this work, open scattering and radiation problems will 
be considered. Consequently, any valid and unique solution must also satisfy the 
Sommerfeld radiation condition [28] which describes the field behavior at infinity 

lim r 

r—f OO 

This simply states that the field is outgoing and of the form e~ jkr jr as r — ► oo. 

Faraday’s Law (2.3) and Ampere- Maxwell’s Law (2.4) are coupled first order 
partial differential equations (PDEs). However, these may be combined to yield a 
single second order PDE in terms of E or H. Specifically, by taking the curl of one 
of these and making use of the other, the following second-order PDE is obtained 


V x 




= 0 


(2.15) 



(2.16) 


These are commonly referred to as the vector wave equations. Solution of the wave 
equation also requires sufficient boundary conditions for uniqueness. 
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Direct solution of Maxwell s equations or the wave equation is only possible for 
a few special conditions. As a consequence, equivalent sources are often employed 
in the solution of practical problems since often such an equivalent problem may 
be readily solved numerically. The surface equivalence principle states that any 
field exterior to a given (possibly fictitious) surface may be exactly represented by 
equivalent surface currents which are allowed to radiate into that external region. 
These equivalent currents are given in terms of a jump discontinuity between the 
total exterior (E\. H\) and interior ( E2 , ^ 2 ) fields 

n x (/fa - H 3 ) = J (2.17) 

nx(£ 2 -£i) = M (2.18) 

The radiated fields due to these equivalent currents are given by the integral expres- 
sions 

E{r) = jj) h(f')x E(f')-V xG{R)dS' - jkZ $h(r') x H(r') -G{R) dS' 

(2.19) 

H(r) = h(r')x H{r')-V xG{R)dS' +jkY ^h{r')x E(r')-G{R)dS' 

( 2 . 20 ) 

where R = |r — r'\, r denotes the observation point and r is the integration point. 
The surface on which the equivalence theorem is applied is denoted by S. In (2.19) 
and (2.20), a dyadic Green’s function is required which satisfies at least the radiation 
condition (2.15). Such a Green’s function is the free-space dyadic Green’s function 

3o = (r +7 T)g<,(b) (2.2i) 

which utilizes the corresponding scalar Green’s function 

-jk 0 R 

G 0 (R) = 


AnR 


( 2 . 22 ) 
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Note that use of (2.21) requires both equivalent currents (2.18). Alternatively, when 
certain bodies such as an infinite metallic plane, cylinder or half plane is present, 
a dyadic Green's function can be chosen which satisfies an additional boundary 
condition such as the Dirchlet condition 


n x Gj = 0 on S 


(2.23) 


so that (2.19) and (2.20) may be written in terms of only one current per field 
expression 


E{r) = n(f') x £(r')-V xG,(fl) dS' 

H{f) = h(r') x H(r')-V xZhiR) dS' 


(2.24) 

(2.25) 


The Neumann boundary condition is also useful 


n x V x G 2 = 0 onS 


(2.26) 


with the resulting field integral expressions in terms of only one current 

E{r) = -jkZ$ s n(r') x H(f‘ ') • G 2 (R) dS' (2.27) 

H(f) = jkY ^h(f')xE{f’)-^ 2 (R)dS' (2.28) 

Such dyadic Green’s functions are commonly referred to as a first kind dyadic Green’s 
function and a second kind dyadic Green’s function, respectively [29]. 

Often, engineering quantities such as the radar cross section (RCS) and the ra- 
diation pattern require evaluation of one of the previous integral expressions as the 
observation point (r) recedes to infinity. For example, the RCS is defined by 


, 2l . 

a = 47t r z lim v — 

r -0O |£«(^| 


(2.29) 


where E“(r) is the scattered field and E l (r) is the incident field. The scattered 
field may be defined by either (2.19), (2.24) or (2.27). Analytical expressions for 
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these fields may sometimes be obtained: however, the required currents are often 
determined via an appropriate numerical method. 

2.2 Numerical Techniques 

Determination of electromagnetic fields within a region typically involve eit her in- 
tegral or partial differential equations. Such formulations may be written in operator 
notation as 


Cu-f = 0 (2.30) 

subject to appropriate boundary or transition conditions 

B (u) = 0 (2.31) 

within the domain (D) and on its boundary (dfi). In these, the operator C is based 
either on an integral representation of the fields such as (2.19)-(2.20) or upon the 
vector wave equation (2.16). The forcing function / is a known quantity while u , 
which is often a field or current, is unknown. Unfortunately, very few analytical 
solutions for (2.30) are available in electromagnetics. One such solution, the fields 
due to a magnetic dipole in the presence of a metallic cylinder, will be used in 
Chapter III to form the appropriate dyadic Green’s function. However, most useful 
electromagnetic scattering and radiation problems cannot be solved using analytical 
methods. Rather, an approximate numerical solution is sought which in some way 
closely resembles the exact solution. Two methods of formulating such an approxi- 
mate solution are: the Rayleigh-Ritz method and the method of weighted residuals. 
Both of these approaches are discussed by Zienkiewicz and Taylor in their classic 
finite element method text [30]. 
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2.2.1 The Rayleigh-Ritz Method 

The Rayleigh-Ritz method seeks a stationary point of a variational functional. 
For operators which are self-adjoint and positive definite, the stationary point of the 
following functional 

F{u) = ^(Cu,u)-{uJ) (2.32) 

is an approximate solution of (2.30). In (2.32), the inner product over the domain ft 
of the two functions is defined as 

(a, b) = [ abdVt (2.33) 

or for vector functions 

(5,6) = ^5 -6 5ft (2-34) 

The choice of this inner product extends the validity of the variational expressions 
to complex fields. The trial function, u, is expanded in terms of N basis functions 

u = ^2 UjWj = {u} T {it;} (2.35) 

i = i 

where wj are the basis functions and Uj are the unknown expansion coefficients. 
In (2.35), column data vectors are denoted with {•} while row data vectors involve 
a transposition Substituting (2.35) into (2.32), the functional may now be 

written 

F(u) = \ {u} T U {u;} C {u;} T 5ft] M - {u} T {w} f dft (2.36) 

This functional is minimized by allowing all partial derivatives with respect to the 
coefficients, {u}, to vanish 

— F(u) = -[/ WiC{w} T dn {u} + l:{u} T [ {w}Cwidn- [ Wif dn 
diii 2 Un J 2 Jn J n 

= 0 


(2.37) 
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which may be written as a matrix system 

MU") = W (->•*<) 

The elements of the matrix [.4] and excitation vector {6} are given by 

A{j = / WiCwj dtt 

Jn 

b t = / wj dn (2.39) 

Jq 

A word of caution, unlike other branches of engineering, no physical significance may 
be attached to the stationary point of the functional (2.32). In mechanical systems 
for example, minimizing this functional represents minimization of the total potential 
energy of the system. However, since electromagnetics involves complex quantities, 
such a statement may not be asserted. Additionally, the variational formulation 
given above requires a symmetric operator C which is not always available in com- 
putational electromagnetics (see for example the FE-ABC formulation in Chapter 
IV). However, the method of weighted residuals does not require a symmetric op- 
erator and is perhaps more familiar to electrical engineers since it does not involve 
variational principles. 

2.2.2 The Method of Weighted Residuals 

The method of weighted residuals attempts to minimize the residual 

n = Cu-f >0 (2.40) 

where 7Z = 0 would imply that u ~ u for all points within Q. Such a solution in 
general cannot be obtained and it is more realistic to find a solution which satisfies 
the minimal residual condition in some average or weighted sense 


J [t t -£ {w} T {u} — $i/] dVt = 0 i = 1,2,3,. ..,N 


(2.41) 
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In general, any testing function (/,-) may be used; however, since these functions 
modify the enforcement of the boundary conditions throughout the domain, the 
choice of testing functions effects the quality of the solution (2.41). One popular 
testing procedure is called collocation or point matching. In this, the testing function 
is a Dirac delta function, t , = 8{x — x,) which implies enforcement of the boundary 
conditions only at discrete points e.g. x t , i = 1,2,3,...,N . Another popular choice 
is termed Galerkin’s procedure. In this, the testing function is identical to the 
expansion function used in (2.35) e.g. t, = w t and the weighted residual equation is 
given by 

yj WiCwj dft, {u} = J Wif dfl (2.42) 

which is identical to the Ritz procedure given above. In Galerkin’s method, the 
boundary conditions are satisfied throughout the computational domain in a least 
squares sense. 

2.2.3 Matrix Assembly 

The construction of the system matrix (2.39) associated with either the Ritz or 
Galerkin’s formulation represents a significant portion of the overall computational 
cost of either the moment method or the finite element method. However, these two 
popular techniques have vastly different cost characteristics which is the primary 
attraction of the finite element method for volume formulations. A general moment 
method matrix requires the computation of all matrix entries which is a 0(N 2 ) 
process where N is the number of unknowns. These matrices are full since each 
unknown couples with each other through the Green’s function. On the other hand, 
FE matrices are sparse due to the locality of the PDE formulation and hence each 
degree of freedom interacts only with its neighbors. Such a matrix assembly process 
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is Q(C A r ) where C is a constant which depends on the type of finite elements used. 
For typical finite elements used in computational electromagnetics, this constant is 
between twenty and forty. Therefore, the matrix build time (and memory require- 
ment) for a FE matrix grows only linearly with the problem size while the build time 
(and memory requirement) for a moment method matrix grows as the square of the 
problem size. Since the number of unknowns increases roughly as the cube of the 
linear dimensions for volume formulations, an FE matrix grows as a cubic function 
of the linear dimension while the corresponding volumetric moment method system 
increases as the sixth power of the same linear dimensions. 



CHAPTER III 


Finite Element-Boundary Integral Method 


The Finite Element (FE) Method is a computational technique which approxi- 
mates a continuous boundary- value problem with a corresponding discrete problem. 
The latter is then solved numerically, thus obtaining a solution which approximates 
in some sense the original continuous problem. The FE method has been used in 
mathematical physics since the 1940’s. It was first coupled with an exact Boundary 
Integral (BI) termination condition in an electromagnetics application by Silvester 
and Hsieh [14] and McDonald and Wexler [15] in 1971 and 1972, respectively. Re- 
cently, there has been renewed interest in the Finite Element-Boundary Integral 
(FE-BI) method of electromagnetics principally due to the work of Jin and Volakis 
[16], [17] and [20] who demonstrated the successful application of the method for 
cavity-backed antennas. Their major contribution was coupling the FE-BI approach 
with the Biconjugate Gradient-Fast Fourier Transform (BiCG-FFT) technique, thus 
allowing high fidelity simulations with a low (9(N) memory and computational de- 
mand. 

This chapter will extend the FE-BI method to cavities which are recessed in a 
metallic circular cylinder of infinite extent. The boundary integral will employ uni- 
form zoning and hence the BiCG-FFT solver may be employed to retain low mem- 
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ory and 0( A log 2 A ) computational burden. We will first present the FF method 
derivation starting with the vector wave equation and then proceeding with the in- 
troduction the boundary integral as an exact mesh truncation condition. An efficient 
computational strategy for the admittance matrix will be presented along with ac- 
curate expressions for interior and exterior source modeling. Finally, the computed 
electric fields will be used to evaluate the far-zone fields and, if applicable, the input 
impedance of patch antennas. 

3.1 FE-BI Formulation 


Consider a cavity recessed in an infinite metallic cylinder, shown in Figure 3.1. 
The cavity walls are assumed to coincide with constant p-, (f>- and z-surfaces and the 
cavity is filled with inhomogeneous material. Interior sources and lumped loads may 
be present as well as surface metallization patterns. 

The FE-BI system is developed directly from the inhomogeneous vector wave 
equation (2.16). To generate a system of equations from (2.16), the method of 
weighted residuals is applied which results in the symmetric inner product of a vector- 
valued weight function and the vector wave equation as described in Chapter II. The 


integro-differential equation is given by 


/ Vx V -Wi{p,<t>,z)pdpd<f>dz 
Jv x fl r [p, 0, ZJ 

~ko / (-r(pA,z)E{p,4>,z)-W l (p,<j),z)pdpd<j)dz = 

J V t 

-/ Vx — 7 ~T } •Wi(p,<t>,z)pdpd<j>dz 
JVi 

—jk o Z 0 J J l {p,(j),z ) • Wi(p, <f>, z)p dp d(j> dz 


where Wi(p,<f>,z) is a subdomain vector-valued weight function to be specified and 
K is the ith volume element resulting from a discretization of the cavity. Given 
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Figure 3.1: Illustration of the cavity geometry situated on a metallic cylinder along 
with the coordinate system used in this work. 
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the impressed sources (J*‘. A/'), the right-hand side of (3.1) is the interior excitation 
function defined by 


/r = - Jh 


V x 


+ jk o Z 0 J'{p,O, z) > • \Y,{p.o. z)pdpdodz{3:2) 


p T {P,<t>,z) 

Upon application of the first vector Green’s theorem [31], we recognize (3.1) as 
the weak form of the wave equation 

/ v x ) • v x mp,M pdf , d<l>ll: 

J \\ P T ( P 1 0 1 ^ ) 

-k\ l e T (p,<f>,z)E{p, 0 ,z) ■ Wi(p,<f>,z)pdpd<f>dz 

Jv, 

-jk 0 Z 0 j>$ n{p , <f>, z) x H(p, <f>, z) ■ Wi{p, <p, z)dS = /•"' (3.3) 


where n(p,<j),z) indicates the outward pointing normal of the element surface asso- 
ciated with the ith unknown, Si is the surface area of that element, and H[p , <£, s) is 
the total magnetic field. It can be shown that the surface integral of (3.3) vanishes 
for all elements which do not border the cavity aperture. Furthermore, their non- 
zero contribution is limited to the portion of their surface which coincides with the 
aperture. 

Unfortunately, (3.3) contains unknown electric and magnetic fields on the surface 
of the cylinder. To eliminate H from (3.3), we utilize the surface equivalence theorem 
introduced in Chapter II (2.28) with a second kind dyadic Green’s function. This 
formula provides a relationship between the electric field and the magnetic field on 
the aperture of the cavity. We expand the total magnetic fields as 


H{r) = H l (r) + H r (r) + H s (f) 
= H cy \r) + H a {f) 


(3.4) 


where H l is the incident field, H r is the reflected field due to the infinite metallic 
cylinder and H s is the scattered field attributed to the cavity. The incident and 
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reflected fields will be specified later in this chapter while the scattered field is rep- 
resented by (2.28). 

Simplifying (3.3) with (2.28) and (3.4). the resulting system is given by 

y ^ £ f f Y. X ^ r ) ' Y. X & ~) 

j 3 JVi \ Pr{p,4>,z) 

-kle r {p, 4>, z)Wj{p, z) ■ Wi(p,<t>,z)^pdpd<t)dz 

+{k 0 a) 2 J2 E : 6 a (j)6 a (i) f f [ll r ,(a, <£, z) • p(a, d>, 2 ) x 

j J S x JSj 

G 2 (a,4),z) x p(a, <f> , z) ■ Wj(a,4> ,z')j d(j> dz'd&dz - f' nt + /f x< (3.5) 
where the unknown electric filed is expanded in terms of subdomain basis functions 


E = J2 E iWj(p,<f>,z) (3.6) 

J 

In (3.6), the subscript indicates the j th unknown and Wj(p,4>, z) is the same edge- 
based expansion function as that used for testing in (3.1); e.g. Galerkin’s procedure. 
The cylinder has radius a and the relative constitutive parameters of the material 
filled cavity are denoted by t r and p r . The function S a (i)8 a (j) is a product of two 
Kronecker functions and it simply indicates that the boundary integral only con- 
tributes when both the test and source unknowns are on the aperture. In addition 
to the interior excitation functional, f- nt , an exterior excitation functional, f - xt , is 
present and will be discussed later in this chapter. The appropriate weight functions 
and the evaluation of the FE matrix entries which are represented by Aij in the 
FE-BI matrix system 



[ 0 ] [01 
[ 0 ] [ 0 ] 


{*?} 






(3.7) 


will be presented next. 
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3.2 Vector Weight Functions 

An important factor in choosing the finite elements for gridding the cavity is the 
element's suitability for satisfying the mathematical requirements of the formulation 
as well as the physical features of the antenna system. Traditional node-based finite 
elements associate the degrees of freedom with the nodal fields but have proven 
unsatisfactory for three-dimensional electromagnetics applications since they do not 
correctly represent the null space of the curl operator and hence spurious modes 
are generated [32]-[34]. In contrast, edge-based elements correctly model the curl 
operator and therefore the electromagnetic fields. In addition, edge-based elements 
avoid explicit specification of the fields at corners where edge conditions may require 
a singularity. Jin and Volakis [20] presented edge-based brick elements which are 
convenient for rectangular-type structures and cavities, but for cavities residing in a 
circular cylinder, shell elements are the natural choice. 

Cylindrical shell elements possess both geometrical fidelity and simplicity for 
cylindrical-rectangular cavities. Figure 3.2 illustrates a typical shell element which 
has eight nodes connected by twelve edges: four edges aligned along each of the 
three orthogonal directions of the cylindrical coordinate system. Each element is 
associated with twelve vector shape functions given by 

Wi 2 (p,<f>,z) = W p (p,<f>,z;-,<l> r ,zt,+), W A z{p,<j>,z) = W p {p, <f>, z\ •,</>*, z t , -) 

W 56 (p, z) = W p (p,<j),z-,-,(t> r ,Zb,-), Ws7(p,(j>,z) = W p (p,<t>,z-,-,(j>i,z b ,+) 

Wi 4 (p, <p, z) = W+(p, 4 > , z-p b , •, z t , + ), W 23 {p, 4>, z) = Wt(p, <j), z; p a , •, z u -) 

W 58 (p,<t>,z) = W^{p,<i),z\pb,-,z b , -), W e7 (p,<f>,z) = W <t> (p,<f>,z\p a ,-,z b ,+) 

W 15 (p, <t>, z) = W z (p,<t>,z;p b ,<t>r,-,+), W 26 (p,4>,z) = W z (p,<}>,z\p a ,<}>r,-,-) 
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1^48 (p. o, ~) — \\ : {p.o. z\ pb.Qi. H’ 37 (p. o, - ) = (3. 


where W tk is associated with the edge which is delimited by local nodes (l.k) as 
shown in Figure 3.2. As seen from (3.8), three fundamental vector weight functions 
are required for the complete representation of the shell element. They are 


Wp(p,<j),z\p, <f>,z,s) 
H<^(p, 0, ^ p, 0 , 2 , 5 ) 
W*(p, 0, 2;p, 0,2,s) 


spb (0- Q)(c - ;) . 

ah p ^ 

ti„ - p )[: - i)i, 

7 ~{p - p){4> - 4>)z 


(3.9) 


where the element parameters (p a , p 6 , 0/, 0r, z b , z t ) are shown in Figure 3.2, t = p 6 -p a , 


a = 0 r - 0/ and h = z t — Zb . Each local edge is distinguished by p, 0, z and s as given 
in (3.8). The ^-term which appears in the definition of the p-directed weight (3.9) 
is essential in satisfying the divergence free requirement, i.e. so that V • Wj = 0 1 . 
Note that as the radius of the cylinder becomes large, the curvature of these elements 
decreases, resulting in weight functions which are functionally similar to the bricks 
presented by Jin and Volakis [20]. Having specified the vector basis functions, we 
may proceed to develop the matrix entries for the system (3.7). 


3.3 Finite Element Matrix 


The inherent locality of a partial differential equation formulation results in a 
highly sparse system, [A]. This FE matrix arises from the discretization of the first 
two integrals of (3.5) and owing to the finite support of the basis functions its entries 

are identically zero unless both the test and source edges are within an element. 

1 Wj(p,4>,z) will only satisfy this requirement within the volume of the element. These weighting 
functions introduce artificial charges on the faces of the element and are not divergenceless at 
element interfaces. This is allowable since these elements do not guarantee normal field continuity 
across the element faces. 
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These matrix entries can be represented as 


a — “ jd)'J L-*, ['* 

-*mj — I st K o t r‘st 

Pr 


,( 2 )tj 


(3.10) 


where we have assumed constant material properties within an element (c r and // r ) 
and the subscripts (i,j) refer to the row and column of the FE matrix, respectively. 
The two auxiliary functions are defined by 


/ ( 1 ) *J 

1 st 


t( 2)0 

l st 


/ V x W a (p, <p , z; Zj,sj) ■ 

Jv, 

V x W t (p,<p,z-,pi, 0 i,zi,si)pdpd(()dz 

f Wsip^^z^pj^j^Zj.Sj) ■ W t (p,(f),z', pi, & ,Zi,$i)p dp d(p dz (3.11) 
Jv, 


where (s,t) G {p,(f>, z} indicate the direction of the source and test edges, respectively. 
Since the fundamental weight vectors (3.9) are aligned along orthogonal directions, 
(3.11) is symmetric with respect to the source and test edges (e.g. l£ ),J = I^' 3 ), 
and thus, only six unique combinations of (s,t) for / ^ need be determined and only 
three such combinations are required for / (2 h They are 


/(i) 

PP 


S s St 

i^hy L 


p\h In 


(-) 

\PaJ J<h 


4>s)(4> - &)d4> + 


“ f^T- 1 ) / {z - Zs){z - Z t )dz 
Zb 


2 \p 2 a 


7 (1) - 

1 p<)> ~ 


Id) 

pz 

l u 


id) _ 


Ss$t 


l 4>z 

iiV 


2p b \n( — ) +pt(l~— )] / (z - z s ){z - z t )dz 
th 2 L \PaJ \ Pa / 1 J Zb 

J_Ap_ I* 

cr J<t>i 

(i “ p °) + \ + ^ ( p “ ” p ’*) + 2 PsPt ~ P “)) + 

^2 [p\ — — 2 1 ( p s + pt) + psPt In ( z ~ z s )(z — z t )dz 

fPb 

/ (p - ps){p - Pt)dp 

J p a 


SjSt 

t 2 j pa 


= ^[«(|W-/4)-‘W. + A) + W. 


5 - Ain u)) + 
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/< 2) = 

Ss*tP b 

In 

pp 

(q/i) 2 

If KJ 

II 

*5 s ^ f Q 

1 

{thy 

.4 



L 

/‘ 2) = 
* zz 

s s s t h 

•i 

(to) 2 

.4 


\(pl-p 2 a )J o (o-o,)(o-o t )rfo] 

f (o - Os){o - <Pt)do I (r - z s )(z - Z,)(h 

J<t>t 


K Pa / 


1 

3 


P&r ~ - 

/ (4> - (p s )((j> - 4> t )d4> 

u, 

where each of the unevaluated integrals are of the form 


(3.12) 


J l (£ - Is) (f “ It) # = 2 (^ 2 “ ^ 2 ) (^ s + I') + 3 (^ 3 ~ ^ 3 ) + 1 *^ “ L) 


(3.13) 


The FE matrix is assembled by evaluating (3.10) for each element combination which 
contains both the i th and j th unknowns. Since the integrals (3.11) are symmetric, 
only the upper triangle or lower triangle of the FE matrix need be computed and 
stored. 

A lumped impedance post may be included in the formulation by adding a term 
to (3.5) and equivalently to (3.10); surface or sub-surface metallization layers may 
also be modeled. Lumped loads are approximated in the FE-BI formulation by a 
filamentary load [35] located at (<f> L , z L ) for radial, ( p L , z L ) for azimuthal and (p L , <j> L ) 
for axial posts, respectively . These posts have length /, cross-sectional area s and 
impedance Zl- Such posts may be represented as 

x t — 4 >l)6(z - zl) , ,. , . 

Z p (p,4>,z) = — — for radial posts 

Zl P 

cx 

Z p (p , 4>, z) = —S{p ~ Pl)&{z - Zl) for azimuthal posts 

Zl 

Z p (p, <j>, z) = ^ — — Pl ^- ~ — — for axial posts (3.14) 

Zl P 
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assuming that each post spans only one element. The contribution to [.4] is given by 

A,j = jk o Z 0 Z p (p,0,z)\\’i(p,0 > z)Wj(p.Q.:)pdpdodz (3.15) 

which may be evaluated in closed form for the three post directions given in (3.14) 
as 








jk.Z.-Li* (^j) 


£-1 


{(d>L - *) [ol - 4>>) ( ~L - 5,-) (z L - £,-)] 


= jkoZoJ- SiSj [(PL - pi) ( PL - P:){z L - =,) (z L - 5 >)] 


Z L 

h_ 

Zl 


jk 0 Z 0 ~SiSj f — j [{p L - pi) {p L - pj) (d> L - 0.) (<t> L - 


(3.16) 


where the superscript indicates the orientation of the post. 

In addition, infinitesimally thin metallization layers may be represented by simply 
fixing a priori the weight coefficients to zero for weights associated with edges which 
are tangential to the metal. This is a consequence of using a total electric field 
formulation. The symmetry and sparsity of the FE system [4] is maintained after 
the addition of these loads while the BI system [Q] remains fully populated and 
symmetric. 


3.4 Boundary Integral Matrix 


The boundary integral provides an exact boundary condition for mesh closure 
and its construction relies on a cylindrical dyadic Green’s function. The entries of 
the boundary integral sub-matrix are 

Gij (fc 0 ct) / / (f>, z, pi, (j>i) Zi, s t ) 

Js, JSj 

p{a, 4 >, z) x G 2 {a, 4>, z) x p(a , <f> , z) 

•W s (a,4> ,z ;pj,4>j,Zj,Sj)d<fi' dz'd<pdz 


(3.17) 
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where the weight functions are given by (3.9) and are evaluated at the surface' /> — a. 
In (3.17). the dyadic Green's function (G 2 ) satisfies both the radiation condition and 
the Neumann boundary condition at p = a. This dyadic Green s function may be 


expressed exactly as in [31] 


G zz (a,frz) = 


r °° / fr 

E / r 

J—oc \ rt a 


W n=— oo 


kX 1 7) 


7 //: <2) (r 


p<bz, 1 r-\ * S" f°° ( n ^ Z ^ 




G^(a,<M) 


1 y f°° 1 [ tfn (2) (7) f V ^n 2) (7) ^(ni-fc.5)^ 

H { n 2) ( 7) V*o«*J //? } (7)J 


(3.18) 


where 7 = k p a and k p = y/k* - k]. However, for large radius cylinders e.g. ka > 3, 
(3.18) is computationally prohibitive. In these cases, which are the main concern of 
this work, it is advantageous to employ an asymptotic expression for G 2 [36]- [39]. 
Such expressions employ a creeping wave series expansion of which only the two 
direct path contributions, as shown in Figure 3.3, will be retained. The formula due 
to Pathak and Wang [36] has proven quite accurate and is given by 

G zz (a , 4>, z) J -t<t e ~ jk ° S { ( Cos20 + 9(1 “ 9)(2 “ W0) ) V{(3) } 

G i ‘ z {a,4>,z) ~ s\n6 cosol (1 - 3<?(1 - q))v{(3) 

Z7T f 

G u {a,j>,z) ~ - J ^qe- jkoS l(sin 2 0 + q{l-q)(2-3sin 2 0j)v((3) 

+q sec 2 $ (u((3) — u(/?))] | (3.19) 

2 

where 3 = ks 3 and q = -j-. In the definition of 3, s is the usual geodesic path 

length (s = ^/(a4>) 2 + z 2 ) and 6 is the direction of the trajectory ( 6 = fan -1 ^ )• 
Depending on which of the two direct paths as shown in Figure 3.3 is used, $ = 4> 
(short path) or $ = 27T - 4> (long path). The soft and hard Fock functions, u((3) 
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Figure 3.3: Geodesic paths on a circular cylinder. 
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and v(3) respectively, are characteristic of on-surface creeping wave interactions and 
have been extensively investigated by Logan [40]. 

Care must be taken in evaluating (3.17) so that the overall storage requirement 
remains 0{ N) and the singular integrals of (3.17) are accurately computed. If uniform 
zoning is used, the resulting sub-matrix ([£]) is block Toeplitz and hence amenable to 
solution using the BiCG-FFT method. For the non-selfcell contributions, mid-point 
integration may be used while a regularization procedure must be employed for the 
self-cell along with a more accurate numerical integration scheme such as four point 
gaussian quadrature. Bird [39] noted that (3.19) recovers the metallic screen Green's 
function when 3 = 0 within the available approximation order. This suggests that 
(3.17) may be regularized by adding and subtracting metallic screen Green’s function 

— - = Wl e~ jkoR , 

2 G,(a,i,z)= ' + -p-“ 2 ^r ! «=l ? - r l (3-20) 

from (3.19). The free-space dyadic Green’s function is given in closed form by (2.21). 
The resulting regularized Green’s function (curvature contribution only) is given by 


G zt (a,hz) ~ - J ^qe-> k ° a {(cos 2 0 + q(l-q)(2-3cos 2 d))[v(0)-l}^ 

G"^(a, <f>, z) ~ ^-qe~ :koS sm 6 cos (1 - 3q(l - q)) [v((3) - 1] j 
G^(a, <j>, z) ~ { sin2 0 + <?(! ~ ~ 3sin 2 d)) b(/3) - 1] 

+q \sec 2 0 (u(/3) - u(/?))] | 


(3.21) 


and since it is not singular it may be numerically calculated. When the regularized 
integrand is used, the contribution of the planar Green’s function 
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p{a.o. z) x G 0 {a.o. z ) x p(a.o . z ) 

■\\’ s ( a. o , z \ pj.Qj, z : . Sj)do dz dodz (3.22) 

is added to (3.21) 

Gij = Gij + G £ (3.23) 


and this is used instead of (3.17) with R = s. Upon use of a common vector identity 
and the divergence theorem [16], we obtain from (3.22) 

( k 0 a) 2 


^ - 27T 


/ / W t (a,<f),z-,pi,(t)i,Zi,Si) • / 

Js x JSj 


• W 5 ( a , 4 > , 2 ; pj , d2 d<^</2 

2tt is Is ^ -2r) x (f > , 2 ? p,*, 2 *, J,)j 


j/cK 


V- [,3(a,</>',z') X W s (a, 4> ,z ] pj,<j>j, zj,Sj)^ - — — d<p dz d(p dz 


(3.24) 


where / = xx + yy + zz. and this form of the boundary integral may be readily 
evaluated even as R vanishes by employing the standard regularization procedure 
(see for example [16] or [41]). This procedure amounts to subdividing the surface 
into smaller patches and using an analytical evaluation of an additional regularized 
integrand. We note that v(/3) — *■ 1 as a — ► oo and hence the regularized integrand 
(3.21) vanishes leaving only the planar contribution (3.24). 

3.5 Excitation: FE-BI 


Two sources of electromagnetic fields are considered: external sources (plane 
waves) for scattering analysis and internal sources (probe-feed) for antenna parameter 
calculations. The use of the exact boundary condition in (3.5) allows the coupling 
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of an exterior excitation field into the cavity while the FE formulation itself readily 
permits modeling of interior sources. In this section we describe the form of the 
source functionals ff xt and f] ni along with their numerical implementation. 

The forcing function, due to exterior sources (/ t er *) is given by 


rexi ■ 

ft = 3 


0 ~° f W x {a, (f) , z ) • p(a , & , r ) x // cy/ (a, <p , r )do dz 

P>ri J S\ 


( 3 . 25 ) 


where p r , is the permittivity of the interior element associated with the i th aperture 
unknown, Wi{p,(f),z) is the corresponding testing weight and H cy represents the 
magnetic field on the cylinder's surface in the absence of the cavity. A plane wave 


E' 


«i e -j ko (k'-r~) 


H' = 


|/5 (sin 7 sin <pi + cos 7 cos 6 , cos 0ij •+- 4> (sin 7 cos (f>i — cos 7 cos #,■ sin 4>ij — 

A * n ?A: 0 a sin 0 t cos 0, cos 0. 

z cos 7 sin u l e J l J 

Y 0 |p (sin 7 cos 6i cos </>,- — cos 7 sin <£,•) — 4> (sin 7 cos 0 , sin </>, -f cos 7 cos — 


z sin 7 sin 6 t 


D jk 0 P sin 6 t cos 0, cos 0* j 


(3.26) 


is assumed to be incident on the cylinder from the direction (#,,<&) where 7 is the 
polarization angle and e l = 0 ‘cos 7 + 0 *sin 7 is the electric field polarization. In 
these, the difference between the observation and incidence angles is denoted by 
4>i — <f) — fa. The total surface field is given by the sum of the incident and the 
corresponding scattered field from the infinite metallic cylinder [42]. Specifically, 


= H i (a,ct>,z) + H a cyl (a,cf>,z) 
= 4>Hf + zHf 


(3.27) 


From traditional modal analysis, we have 

o jk 0 cos 6 t z 00 


H?M,z) = -2 Y 0 


cos 7 


nk 0 a sin 0 t n fr ' 00 [m 2) (k o asrn0 t ) 
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n sin -v cos 0 X 

k Q o sin 6 X Hn 2 \k 0 a sin 0 ,) 

H c z yl (a, <Z>, z) = j2Y 0 S -^e jk ° cose ' : £ 

7T k n Q 


r jn(^ + C-0, ) 


e Jn(f +0-0, ] 


Hn 2 \k 0 a s\n6 t )\ 


(3.28) 


These expressions may be approximated by retaining only a few terms of the series 

if kodsmO, is small. However, as this parameter becomes large (e.g. for large a and 

0 , — > 90°), (3.28) may be replaced with equivalent asymptotic representations similar 

to those considered earlier. Utilizing Watson’s transformation and Fock theory [42] 

in connection with (3.28), we find the following direct path creeping wave terms 

2 

Hf ~ —Y 0 sin 7 sin ^e jfc ° cosfl,z ^ e -jk 0 as\n6,<t p jg(°)( m <|> p )] * 


P-1 


Hf ~ j2Y a cos 7 m . - e ’ k ° cose ‘ z Y e~ 3k ° asin8 '*r 
9 k 0 asin&i L J 


p= i 


-Y 0 sin 7 cos 0 t e jkoCOs9 ' z £(-l 

P - 1 


S (0, (m$ p ) 


-J r 

k 0 a sin 0 , 


(3.29) 


in which $i = — ( 4 > — <f>,), $> 2 = ( 4 > ~ <t>i) — f , m = ] 3 , and the asterisk 

indicates complex conjugation. The appropriate Fock functions are [40] 2 


(3.30) 


9 "’ (f) “ ^Irw l (t) dt 

f m (() = -4= / —rr.dt 

where Wi(t) denotes the Airy function of the first kind, te^) is its derivative with 
respect to the argument and the integration contour (T) is given by Logan [40]. 

The asymptotic formulas (3.29) are quite accurate except when 4> K 4>i- In this 
region, Goriainov’s [43] expressions 

Hf ~ —Y 0 sin a sin 0 I e jA '° cos ^ z Ie -jA;<>asinei * 1 [s ( 0 ) ( m$ i)]* 


2 Logan [40] uses the e time dependency in the definition of these functions requiring the 
complex conjugation in (3.29) 



36 


+ C jkoas ' n6 ' cos( ^- 0l) [<S*(— 77i C0S(0 - O,))]* | 
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H <t> ~ 
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Ar 0 o sin 0 X 
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k Q a sin 
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G ,(,) (e) = g (,) (ty" 

F {1 \ 0 = / (l) (0e*'£ 


(3.31) 


(3.32) 


have been found to be more accurate and can be used instead of (3.29). These 
surface field expressions may now be used to calculate the entries of the column 
vector {fi Xt } via a numerical evaluation of (3.25). In particular, the modal series 
(3.28) is used when k 0 as’m6i < 10 and either (3.29) or (3.31) when k o as\n0i > 10 as 
appropriate. Having now specified the external source, we present the more compact 
internal source used for this work. 

Conformal antenna patches are typically fed by a microstrip line printed along 
with the radiator on the surface of the substrate. The microstrip lines are in turn fed 
by a coaxial probe which originates behind the cavity as shown in Figure 3.4. The 
patch may also be fed directly by the coax feed or through some form of aperture 
coupling. Nevertheless, the principal excitation for the system is given by (3.2). The 
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impressed excitation current for an infinitesimally thin probe is given as 


J' nt (p. <?,-) 

J int (p,4>,z) 


Ph 


8(o - Of)6(z - zj) 

P 


oioHp - - z j) 


8(p - pj) 8(4>- d>f) 


for radial probe 

for azimuthal probe 

for axial probe (3.33) 


which results in an excitation function (3.2) 





J ko Zo I o 


S ' Pa In 


ah 



(<t>j - 4>i) ( z j - ^)] 


for radial probe 


-jkoZ 0 I 0 ^ \{P} Pi) i z f ^t)] 

- jkoZ 0 I 0 

ta 


\(Pf - Pi) {<t>f ~ <?>•)] 


for azimuthal probe 
for axial probe (3.34) 


With both types of excitation and the FE-BI matrix now specified, the BiCG method 
may be used to determine the unknown electric fields within the cavity. 


3.6 System Solution 


The FE-BI system (3.7) may be solved using one of the popular direct meth- 
ods such as Gauss-Jordan elimination, Gaussian elimination and LU decomposition 
[44]. Alternatively an iterative method such as the Gauss-Seidel method [44], the 
conjugate gradient method (CG) [44] or the biconjugate gradient (BiCG) method 
[44]- [46] may be used. The symmetric form of the BiCG solver has been chosen since 
it requires much less memory than a direct method and it can be implemented in a 
manner which is computationally efficient; i.e. utilizing only one matrix-vector prod- 
uct per iteration. Although use of an iterative method such as the BiCG method 
can require more wall clock time than a direct solver when multiple right-hand sides 
are considered, in-core memory is deemed the most critical and expensive resource. 
The BiCG algorithm is given in Appendix B and the efficient FFT-based calculation 
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of the boundary integral matrix-vector product is discussed in Appendix C. In this 
section, we will present only details of this iterative solver specific to this application. 

The BiCG algorithm requires one matrix- vector product per iteration (see the 
algorithm in Appendix B). This operation represents the bulk of the computational 
demand of the method and requires 0(N 2 ) complex operations per iteration for fully- 
populated matrices. The matrix-vector product is carried out by executing the sum 

N 

y[n] = [A]{x) = £ A[n,n']a:[n'] n = 1,2,3, ..., N (3.35) 

n'=l 

and if the matrix is sparse, a storage scheme such as the Compressed Sparse Row 
(CSR) format may be used to reduce both the memory demand and computational 
load. The FE matrix [A] in (3.7) is such a sparse matrix. The CSR scheme retains 
only the non-zero entries of the matrix in one long data vector along with another 
data vector, known as the offset vector, which contains the number of non-zero entries 
per row of the matrix. Using the CSR scheme, (3.35) can be rewritten as 

y[n) = [A] {x} = ^2 A |e(n,n')J x[n] n = 1, 2, 3 , N (3.36) 

n'=l 

where r[n] is the number of non-zero entries per row of the matrix and e(n,n ) 
indicates which entry of the long data vector is associated with the matrix entry 
A[n, n']. Typically, the rows of the matrix are stored consecutively which results in a 
counter being used for the mapping function e(n,n). Although additional memory 
saving is possible since the FE matrix is symmetric (thus only the upper triangle of 
the matrix need be stored), experience has shown that use of a symmetric matrix- 
vector product leads to a severe performance degradation on vector computers (such 
as a CRAY) due to the resulting short vector lengths. Therefore, the entire sparse 
FEM matrix is stored and used in the product so that the code’s performance is 
maximized when executed on vector architectures. 


40 


The boundary integral matrix-vector product involves the fully populated matrix. 
[Q]. If uniform surface elements are used in the discretization, this matrix- vector 
product may be expressed as a truncated discrete linear convolution and is thus 
amenable to efficient calculation using the fast Fourier transform (FFT). Although 
uniform zoning imposes restrictions on the geometries which can be analyzed by this 
FE-BI technique, the resulting memory and computational efficiencies have proven 
to be well worth this sacrifice. 

The boundary integral product is implemented as described in Appendix C with 
the following problem specific exception. The cross-term arrays do not possess the 
property: G^ z [m — m\ n — n] = G^m — m,n — n\ and hence the periodic replica- 
tion rule used by Jin and Volakis [47] cannot be used here. Replication refers to the 
padding process necessary in the computation of a linear convolution via a circular 
convolution. A different replication rule was developed using traditional techniques 
to perform the discrete convolution which is described in Appendix C. A general dis- 
cussion of replication is provided in that appendix. Cylindrical-rectangular cavities 
are bounded by metal on all four lateral sides while wraparound cavities encircle the 
cylinder as one unbroken groove. These two types of cavities are shown in Figure 
5.7. The cross-term replication rule must be changed for wraparound arrays since 
the cavity discretization is fundamentally different than is the case for cylindrical- 
rectangular cavities. We note that these replication rules are not unique and are 
implementation dependent. In particular, since the columns of the discretization 
wraparound on a circular cylinder forming closed loops, the replication rule must 
account for this periodicity to reduce the computational burden. 
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3.7 Scattering and Radiation 

Once the cavity aperture and volume electric fields have been determined for 
either an external excitation (scattering) or an internal excitation (radiation), sev- 
eral engineering quantities may be calculated. The aperture fields may be used to 
determine the radar cross section (RCS) for scattering or the radiated fields for an- 
tennas. This entails the integration of the aperture fields with an appropriate Green's 
function. In addition, the input impedance may be calculated by using the interior 
cavity fields. In this section we will present the relevant formulae for calculating the 
far-zone radiated/scattered fields and the input impedance from the electric fields. 

3.7.1 Far- field Evaluation 

Two of the most important applications of the presented formulation deal with 
the calculation of the cavity’s RCS and the radiation pattern due to sources placed 
within or on the aperture of the cavity. 

We begin with the integral representation of the scattered magnetic field in terms 
of the aperture fields. We have, 

H s ( r ,6,<j>) = jY 0 k 0 a J^G 2 {r,0,faa,<f> ,z') ■ 

[/5(a, <f> , z ) x E(a , <f> , 2 /)] d4> dz (3.37) 

with (r indicating the observation point in spherical coordinates. When the 
observation point is very far from the cylinder, the dyadic Green’s function in (3.37) 
can be replaced by its far-zone representation 

-jk 0 r 

<5 2 (r, 9, fa a , <f>\ z ) ~ \G h 0$ + G ez 0z + G*+H>' } (3.38) 

rC 0 r 

where the unprimed unit vectors are functions of the observation position and primed 
ones are functions of the integration point in (3.37). The components of this far-zone 
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Green's function are determined by a mode matching procedure giving 
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(3.39) 


As one might expect, these series converge rather slowly for large A: 0 a sin They 
must therefore be recast in another form by employing Watson’s transformation and 
Fock theory as described before. We have, 
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p=l 


where the appropriate Fock functions are given by (3.30). As expected, when <f> « 
the formulas attributed to Goriainov [43] 
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are more useful. The far-zone scattered or radiated field can be computed numerically 
by using either (3. 39), (3. 40), or (3.41) in (3.37). 
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For the scattering problem, the RCS is most often the quantity of interest which is 
given by (2.29). Alternatively, the antenna gain may be computed from the far-zone 
fields as 


GdB{Qi 4 >) = lOlogjo 


47T 


fc) ; 


|£ r (M)P 


+ 101og 10 


r a 
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cm 


iZ 0 R xn \ 


(3.42) 


where A cm is the wavelength in centimeters, R tn is the input resistance which is given 
in the next section and E T is the radiated electric field as r — ► oo. In the far-zone, 
the electric field is related to the magnetic field by 


Ej, = —Z 0 H$ 

E e = Z 0 H * (3.43) 

3.7.2 Input Impedance 

In addition to the antenna pattern and gain, designers are concerned with the 
input impedance of an antenna for feedline matching purposes. The FE approach 
allows the calculation of the input impedance of the radiating structure in a rather 
elegant manner. The input impedance is comprised of two contributions [48] 


Zp + Zp> 


(3.44) 


where the first term is the probe’s self-impedance which is the impedance of the 
probe in the absence of the patch and the second term is the contribution of the 
patch current to the total input impedance. The self-impedance of a radial probe is 
typically approximated as [49] 


Zp = 60 


(k 0 l) 2 y/r r + jkj sin 
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-4 


r ‘ + 4/ 2 


2 / 


(3.45) 


where t r is the relative permittivity of the substrate, r c is the radius of the probe 
and / is the length of the probe. For very thin substrates and thin probe wires, this 
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contribution is negligible. Ignoring the probe-feed's self impedance, we have [ IS] 

Z\ n = -yj / E(p.O.z) - J\ nt {p.O.:)pdpdod: (3.46) 

where the impressed current is given by (3.33), Vy refers to the volume elements 
containing the probe-feed, the electric field is the interior field associated with the 
feed edge and I a is the constant current impressed on the probe. Utilizing (3.9) and 
(3.33) into (3.46) yields 


zL = 

l St ^ ln ( _ ) (^ S ^*) 

lo ah \PaJ lK 7 J 

for radial probes 

Z\n = 

j Si th ^P s 

for azimuthal probes 

Z\n = 

/ ^ Si ta [( Ps ^)1 

for axial probes 


which must be summed over the four parallel edges of the element which contains 
the feed. In (3.47), the probe location is indicated by ((f> s ,z s ), ( p 3 ,z s ) or (p s ,<f> s ) for 
radial, azimuthal or axial probes, respectively. 



CHAPTER IV 


Finite Element- Absorbing Boundary Condition 

Method 


In a previous chapter, the Finite Element-Boundary Integral (FE-BI) Method 
suitable for cavity-backed antenna analysis was introduced for cavities which are 
recessed in an infinite, metallic cylinder. The interior fields of the cavity are deter- 
mined via the FE method while the aperture fields are determined using a boundary 
integral in conjunction with the FE method. The boundary integral provided an 
exact mesh closure condition; however, use of such a global condition partially de- 
stroys the sparsity of the FE-BI system. In fact, as the surface area of the aperture 
increases, the poor scaling property of the boundary integral dominates the compu- 
tational and memory demand of this methodology. In addition, the FE-BI method is 
ill-suited for use with coated conformal antennas or protruding antennas. Therefore, 
it is advantageous to seek an alternative mesh closure condition which preserves the 
inherent sparsity of the FE system. 

One approximate closure condition utilizes an absorbing boundary condition 
(ABC). This condition provides an approximate relationship between the tangential 
electric and magnetic fields on the boundary of the computational domain. Wilcox 
[50] introduced a spherical wave expansion of these fields which is the typical ap- 
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proach taken for deriving ABCs. A major disadvantage of two popular ABCs for 
the vector wave equation, which were proposed by Peterson [22] and Webb and Kan- 
nelopolous [23], is the requirement that the mesh boundary be spherical. For most 
geometries, including cylindrical-rectangular cavities, the use of a spherical mesh 
closure results in an excessive number of unknowns. Consequently, conformal ABCs 
which minimize the number of unknowns expended in the region between the geom- 
etry and the outer boundary have been actively sought. One conformal ABC which 
has proven remarkably accurate for scattering applications was recently proposed by 
Chatterjee and Volakis [24]. This ABC uses a second-order expansion of the local 
fields on the boundary to relate the electric and magnetic fields and it explicitly 
includes the curvature of the boundary in its formulation. 

In this chapter, the FE method will be married with these new ABCs to form a 
FE-ABC method suitable for cavity-backed antennas which are recessed in a circular 
cylinder. The complete formulation requires development of the FE equations along 
with the specification of the ABC for a cylindrical-rectangular box boundary. An 
important aspect of this development is the use of a domain decomposition which 
employs a scattered field formulation in the exterior region while utilizing a total 
field formulation within the cavity. Such an approach minimizes the number of 
unknowns within the cavity while removing the incident and reflected fields from 
the exterior region unknowns. The ABC presented by Chatterjee and Volakis [24] 
absorbs only the scattered field in the exterior region (see Figure 4.1) and due to 
known error propagation effects associated with a total field formulation, a scattered 
field formulation is preferred in the exterior region. For radiation analysis, such a 
decomposition is unnecessary since the radiated field is identical to the total field. 
Previously, the FE-ABC method had only been used for scattering analysis of finite 
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bodies (see for example Chatterjee ei al. [21]). A primary goal of this chapter is to 
present the first application of this method to infinite bodies as well as the first use 
of ABCs for antenna parameter calculations. We begin with the formulation of the 
F E-ABC equations. 

4.1 Formulation 

Consider the computational domain shown in Figure 4.1. There are two volume 



Figure 4.1: Typical coated cavity-backed patch antenna with ABC mesh termination. 

regions: an exterior region, V 1 , which includes any radome overlay and an interior 
region, V 11 . Both regions may be inhomogeneous and are separated by the aperture 
surface, S ap , and the surface metallization surface, S 3m both of which lie on the 
surface of the metallic cylinder ( p = a). Thus, the exterior region is defined by p > a 
while the interior region has p < a. The computational domain is bounded by the 
union of the metallic surface, S metal = S sm + S'™ where S cm is the metallic walls of 
the cavity and the ABC surface, S abc . 

Within the computational volume, the total electric fields may be written as 


E(r) = E\r) + E cy \r) f e V 1 



4S 


= E u {r) re V" H.M 

where E cyl (f) = E'(r ) + E T {r) as before. The total magnetic fields are likewise 
written 

H(r) = H l {f) + H cy! {r) f e V 1 

= H U {r) reV 11 (4.2) 

where H cyl {r) = H'{r) + H T {r). The boundary conditions are readily written in 
terms of the electric and magnetic fields. Within the cavity, the tangential electric 
field vanishes on the metallic walls 

hxE n {r) = 0 reS cm (4.3) 

while on the aperture, the total tangential fields are continuous 

n x E\r) = n x E n (r) r € S ap 

n x H\f) = fix H n (f) -fix H cyl (r) f € S ap (4.4) 

On metallic surfaces, all tangential electric fields vanish, i.e. 

hxE cyl {f) = hxE^f) = n x E n (f) = 0 fer (4.5) 

while ii x E cyl (f ) also vanishes over the aperture 

h x E cy \r) = 0 r € S ap (4.6) 

since it contains both the incident and reflected fields. Thus, the only non-zero 
electric fields on the surface of the metallic cylinder correspond to the unknown 
fields within each region which are continuous across the surface aperture as implied 
by (4.4). 
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The FE equations may be developed by considering the inhomogeneous vector 
wave equation (2.16). Employing the method of weighted residuals and Green's iirst 
vector identity, the weak form of the vector wave equation (3.3) is obtained where the 
interior source functional, f \ nt , is once again given by (3.2). A domain decomposition 
is accomplished by substituting the total field relationships (4.1) and (4.2) into (3.3) 
producing the FE-ABC equation 
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However, due to the boundary conditions at the surface of the cylinder, (4.4) and 
(4.5), the two surface integrals over the aperture cancel one another and with some 
further manipulation, (4.7) can be written as 
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(4.8) 


where we have used Ampere- Maxwell’s Law (2.4) to express the magnetic field on 
the ABC surface in terms of the curl of the electric field. The second term on the 
right-hand side of (4.8) serves the same role as the exterior source functional, / t ex< , 
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given previously for the FE-BI formulation in (3.25): however, as will he shown, this 
term will vanish when combined with the succeeding volume contribution. The last 
term on the right-hand side is a source functional which is present due to the use 
of a scattered field formulation in the exterior region, V 7 . As written in (4.8), this 
final term must be computed over the entire exterior region. However, by applying 
Green’s first vector identity, this final term may be written as 



The quantity within the brackets in the volume integral of (4.9) is recognized as the 
vector wave equation and since E cyl already satisfies the free-space wave equation, 
this contribution is non-zero only for material ( c r ^ 1 and/or \i r ^ 1) regions. The 
surface integral is the boundary term which is non-zero only over the ABC surface, 
the aperture surface and the interface between two materials in the exterior region 
which have different permeabilities. Thus, (4.9) may be rewritten as 
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M J sr (ft X ff*') -W.dS+j^ J s;r (, 5 x /?“") ■ W, dS + 

- Lj / n • \Wi x V x E cyl 1 dS (4.10) 

Jiri ftre J Jsj d 


where /i ri is associated with the interior material layer and fi re is the corresponding 
quantity for the exterior material layer. The third surface integral need only be 
computed for each dissimilar material interface in the exterior region. Utilizing the 
source-free vector wave equation and assuming that material properties are constant 
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within an element, (4.10) may be included into (4.8) to get the final FE-ABC equation 


X, 

Iv n 


V xE' YxW, 
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J s id 



dS 

(4.11) 


The source contribution associated with the ABC surface has been canceled by the 
addition of (4.10) and the remaining terms of the right-hand side of (4.11) are consis- 
tent with scattered field formulations. The aperture contribution is analogous to the 
exterior source functional used for the FE-BI formulation with the exception that 
the permeability is now associated with the finite element on the exterior side of the 
aperture. Implementation of (4.11) will give superior results to (4.8) due to numeri- 
cal error propagation which is inevitable when the volume source term is evaluated 
throughout the exterior region as in (4.8). In addition, since (4.8) need be computed 
over the entire exterior region rather than over the presumably smaller material re- 
gion, it would require substantially more computational effort than is required for 
(4.11). This later set of FE-ABC equations may be written as a linear system of 
equations 



e—abc 


{£*"} 



(4.12) 


where the FE matrix [Af e -abc] may be written as a sum of the FE matrix used in 
the FE-BI formulation (3.10) and a second term attributed to the ABC surface 



[*^/e— afcc] — [«4] T 


(4.13) 
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The new FE-ABC equation (4.11) is comparable to (3.3) except that the latter 
utilizes a total field formulation throughout the computational domain. However, 
previously we utilized an integral expression for the total magnetic field across the 
aperture which resulted in the FE-BI equation (3.5). Such an integral expression pro- 
vides an exact relationship between the total tangential electric and magnetic fields 
over the aperture surface w'hich also formed the computational domain boundary. 
Alternatively, we may employ an approximate relationship between these two fields 
with the goal of retaining the sparsity of the resulting linear system. Additionally, 
as shown in (4.11), this FE-ABC formulation may be used for coated as well as un- 
coated geometries. In the next section, we will develop an approximate relationship 
suitable for mesh closure. 

4.2 Conformal ABCs 

Traditional three-dimensional vector ABCs [22, 23] require a spherical outer 
boundary which results in an excessive number of unknowns. New conformal ABCs 
have recently been proposed by Chatterjee and Volakis [24] which have an outer 
boundary that follows the contour of the enclosed geometry resulting in a minimal 
number of unknowns. In this section, the specific expressions required by this new 
ABC for a cylindrical-rectangular box boundary will be derived. A definition of ABC 
order will be given and subsequently the first- and second-order ABC expressions will 
be presented. 

For the purposes of discussion, we define a secondary field as the field which 
is a consequence of equivalent currents that are supported by some primary source 
which is either external or internal to the computational domain. Thus for scattering 
problems, the scattered field is the secondary field while the incident and reflected 
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fields are considered primary fields. Likewise, for a radiation problem, the radiated 
field is the secondary field whereas the source field due to an impressed current is 
the primary field. In (4.11), we recognize that an ABC must supply a relationship 
between the tangential components of magnetic and electric secondary fields on the 
absorbing boundary, S abc . 

The secondary field may be expressed as a Wilcox expansion 


E s {n,t\,t 2 ) — 


D -jk 0 n 


E 


Ep{t i, t 2 ) 


47 Ty/u P^oo ““ U p 


(4.14) 


where u = \/R] # 2 , = Pi + n and pi is a principal radius of curvature. In this 

form, the curvature of the non-spherical wavefront is explicitly used. The point of 
observation is given in Dupin coordinates [31] as 


x = nh + x 0 (ti,t 2 ) 


(4.15) 


where h is the unit normal and x 0 (ti,t 2 ) denotes the surface of the reference phase 
front and therefore, t\ and t 2 denote tangential coordinates on that surface. Absorb- 
ing boundary conditions annihilate outward propagating waves up to a certain order. 
A zeroth-order (P = 0) ABC represents the usual Sommerfeld radiation condition 
(2.15). A first-order ABC (P = 1) annihilates all fields with up to a u -1 dependency 
while all higher order fields are reflected back into the computational domain. For 
a cylindrical surface, it = v /p, thus the zeroth-order ABC is simply the geometrical 
optics spread factor while the first order ABC annihilates fields up to 0(p~ l ). Ev- 
idently, as the ABC order increases, the reflected fields have an increasingly higher 
attenuation factor and hence the boundary may be placed closer to the geometry 
without inducing erroneous reflections. 

In the following two sections, we will present the first- and second-order confor- 
mal ABCs attributed to Chatterjee and Volakis[24]. In particular, the appropriate 



expressions for a cylindrical-rectangular box boundary will be given. The matrix en- 
tries which enforce these ABCs will be presented assuming the use of the edge- based 
weight functions (3.9). 

4.2.1 First-order ABC 


The first-order ABC is similar to the impedance boundary condition attributed to 
Rytov [51]. Specification of the ABC requires the relationship between the tangential 
components of the electric field and its curl. For the first-order ABC. we have 


n x V x E s = \jk Q + K m — ?p] 


(4.16) 


— * A A ___ A A. A A 

where E l — t\E T t 2 El i ^ n m — r/ — 4 - /V2^2^2 snd /ci t 2 are the principal 

curvatures of the surface. These curvatures are given by [31] 


1 dh i 
hi dn 


(4.17) 


where ^ denotes a normal derivative and /i, are the appropriate metric coefficients. 

The ABC formulae are going to differ for each surface of the cylindrical-rectangular 
box. There is a singly curved surface which has h — p with p being constant and a 
total of four flat surfaces having either h = with </> being constant or n = ±z with 
z being constant. Since the curvature and tangential edge orientations are different 
for these three cases, each surface will be examined separately. 

For singly curved surfaces, the normal unit vector is radially directed and the 
tangential edges which form the surface are cf>- and i-directed. Accordingly, (4.16) 
may be written as 


p x V x E s = 


jK + 2 \ 




ik °-h 


Eli 


(4.18) 


since = <f>, t 2 = £, Ki = — - and k 2 = 0. The matrix entry for the ABC surface is 
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given by 


,4 afc<r = J [n x V x Wj] • U , dS 


(- 1 . 19 ) 


as extracted from (4.11). For the singly curved surface, dS = pdcpdz and thus 
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( 4 . 20 ) 


where the subscripts denote the orientation of the test and source edges, respectively, 
and the edge-based weights (3.9) are used for IV, and VF, in (4.19). 

For the flat surfaces which have h = ±0, the ABC (4.16) is given by 


0 x V x E* = jk a [pE s p + zE* 


(4.21) 


since the tangential edges are oriented in the p- and £-directions. The corresponding 
matrix entries (4.19) are 



For the flat surfaces which have n = ±z, the ABC (4.16) is given by 


zxVxE* = jk 0 [pE s p + 4>e; 


(4.23) 


since the tangential edges are oriented in the p- and ^-directions. The corresponding 
matrix entries (4.19) are 





[3 ~~ $*) — 2 (^’ _ $) + 

pi) - | (pi + Pj) ( Pb - P 3 a) + \piPi (pb - Pa)] 


(4.24) 
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These ABC matrix contributions may be assembled in the usual fashion and 
added to the FE matrix (3.10) previously used with the FE-B1 formulation. Since 
the preceding ABC has only first-order accuracy, the mesh closure must often be 
placed far from the cavity aperture which results in a large number of unknowns. To 
reduce the computational burden of the method, second-order conformal ABCs need 
to be explored as presented in the next section. 

4-2.2 Second-order ABC 

The conformal second-order ABC is necessarily more complex than the first-order 
ABC. The tangential component of the curl may be approximated as 

n x V x £ s = g-£ t 4 + /l-Vx [h(n-Vx£ s )] + f-V t (h-£ s ) (4.25) 

where denotes the tangential surface gradient operator [31]. Unfortunately, use 
of (4.25) would result in an asymmetric system [,4 a6c j due to the last term which 
possesses only one differential operator. An asymmetric system requires an itera- 
tive solver which utilizes two vector-matrix products such as the conjugate gradient 
squared (CSG) solver presented in Appendix B which is based upon the work of 
Vorst [52]. A symmetric system, as seen previously in Chapter III, requires only 
one matrix-vector product if the BiCG solver is used. Additionally, for symmetric 
systems, only the upper or lower triangle of the matrix need be computed and stored. 

Chatterjee and Volakis [24] realized that the gradient in (4.25) may be approxi- 
mated by 

V t (V-E 3 t ) = jk 0 V t (h-E°)+0(n~ 5 ) (4.26) 

With both a gradient and a divergence operator present, one operator can be trans- 
ferred to the test vector while the other may remain with the source vector. Hence, 
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the resulting matrix may be symmetric since both the test and source fields are 
differentiated. With (4.26). (4.25) may be written 


fixVxE 1 = ^-E’ + 3-V x [n(n-Vx E s )] + vV,(V-£ ( s ) (4.27) 


For the basis vectors (3.9) considered in this thesis, V • Ef is always zero on S abc and 
hence the third term of (4.27) will not contribute to this form of the ABC. For surfaces 
with a common constant curvature for both tangential directions on a surface, this 
new vector ABC (4.27) will lead to a symmetric FE system j A abc as shown by 
Chatterjee [25]. However, if the principal curvatures on a surface are unequal, the 
system will be asymmetric. For either (4.25) or (4.27), the three coefficient dyads 
are given by 
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UU 


(4.28) 


where k „ = K1K2, Ak = ki — k 2 and D = jko + 5 K m — ^ SL . In the case of (4.27), 7 
must be divided by jk Q due to (4.26). 

As was the case for the first-order ABC, it is advantageous to consider the second- 
order ABC for singly curved and flat surfaces separately. For a singly curved surface, 
the unit normal direction and curvature parameters are fi = p, n 1 = — ^ and k 2 = 0 
as before. After some manipulation, we find that (4.28) becomes 
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(4.29) 
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where D — jk Q — Note that 3 is not symmetric unless /?—►:> c 

For the second-order ABC, it is advantageous to segment the matrix entry (-1. 19) 
into three parts. After some vector manipulation [25], these contributions are given 
by 


j(l)abc 
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dS 
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(4.30) 


where /( 3 ) abc is only used for the original, asymmetric ABC (4.25). On the curved 
surface, we find 
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(4.31) 


where the first subscript denotes the test edge direction while the second subscript 
indicates the source edge direction. Note that the cross terms of the second contri- 
bution has I^J abc 7 ^ I^ abc which results in an asymmetric FE system regardless of 
which form of the ABC, (4.25) or (4.27), is used. This asymmetry is a direct conse- 
quence of the different principal radii of curvature for a cylindrical surface. However, 
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these terms are asymptotically identical as the radius of the ABC surface becomes 
large since the surface will then be approximately planar. A symmetric ABC may 
be obtained by dividing both the numerator and the denominator of (4.31) by the 
ABC radius, pt which results in 


r(2)abc _ 

l <t>z 


j(2)abc _ 
1 z<t> 


(4.32) 


To determine whether the symmetric evaluations given above may be used in place 
of the asymmetric terms, the RCS was computed using all three forms of the ABC 
for a 0.6A o x 0.6A o cavity recessed in a 0.2A o metallic cylinder. The cavity if filled 
with a material which has t T = 2. The ABC surface had a radius of 0.5A o and the 
lateral walls were place 0.3A o away from the cavity sides. Thus, the ABC surface is 
placed at least 0.3A o away from the cavity aperture which is shown to be sufficient in 
Chapter V. The computed RCS for all three ABCs were within 0.01 dB of each other 
and 0.5 dB of the value computed by the FE-BI method. This comparison shows 
that the symmetric form of the ABC may be used in place of the asymmetric ABC 
in order to reduce the computational and storage demand of this FE-ABC method 
even for a highly curved surface. 

For the flat surfaces (n = ±0 or n = ±z ), we find that the first contribution to 
the ABC entry, l( 1 ) abc ^ is identical to the first-order case, (4.21) or (4.23), since 


a = jk 0 \titi +t 2 t 2 


(4.33) 


where i\ = z, i 2 = p if n = ±4> or t\ = p, t 2 = <j> if n = ±z. However, for the 
second-order ABC, there are additional contributions due to 


/3 = -rj— [Mi + M 2 

J K'O 

7 = jKP 


(4.34) 
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However, if n = ±z, these additional contributions are given by 

r(2)afec _ ( Pb\ _ i 

pp ~ j2k 0 a [\ pj 

k 0 t 

k Q t 

nr = ^(?)[*wv.)-2(ft + «x + ««i"(a)] 

/»?*'" = SiSj (^7) [j (ft ~ft)~\ (f + fi) (ft ~ ft) + ifi° 
i£’ tc = + < 4 - 36 ) 

Although it appears that I^ abc / /^ a6c for n = ±z, in fact these terms will be 
identical when evaluated on the ABC surface since /5, = p 3 = p a . 

We have now fully specified the admittance matrix for both first- and second- 
order ABC mesh closure conditions. When coupled with the previously given FE 
matrix, the resulting system retains (9(N) sparsity regardless of the problem size 



or composition. Furthermore, coated geometries may be readily modeled with this 
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formulation since the FE method permits inhomogeneous materials. However, the 
source terms given previously as the right-hand side of (4.11) must still be specified. 

4.3 Excitation: FE-ABC 

In the previous section, vector ABCs were presented along with the associated 
matrix entries required to implement the FE-ABC formulation (4.11). We can now 
generate the FE matrix, [Af e - a bc], which is fully sparse. In order to find the scattered 
or radiated field, we need to develop the excitation function (or right-hand side) of 
(4.11). The first term, is the interior current probe source given previously by 
(3.2). The second term is analogous to the external source functional, /, eit , used for 
the FE-BI formulation(3.25); however, this term now contains the permeability of 
the material exterior to the cavity aperture rather than that of the cavity fill. The 
final term is a consequence of the scattered field formulation in the exterior region. 
It vanishes for free-space since the incident and reflected fields satisfy the vector 
wave equation (2.16). For coated geometries, this term will not vanish and must be 
retained for that case as discussed later. 

4.3.1 Aperture Term 

Previously, the FE-BI formulation required the external excitation function given 
by (3.25). In this, the magnetic field attributed to the incident and reflected fields 
is evaluated on the mesh boundary. Since the mesh boundary coincided with the 
cylinder surface, considerable simplification was possible in generating these magnetic 
fields (see (3. 28), (3. 29) and (3.31)). However, for an FE-ABC formulation, the mesh 
boundary does not correspond to the cylinder surface although as shown above, a 
source term may be attributed to the aperture. Therefore, in the absence of any 
radome covering, the exterior source term is similar to the one used for the FE-BI 
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formulation (3.25). 

The forcing function due to exterior sources is given by 

fext(abc) _ -^aZo f ^ ( a i0 '. ) . f, ( a , q . ) x H cyl ( a . 6 . c' ) do' dr ' (4.37) 

‘ Pre 

where the integration surface, S„ is associated with the i th unknown which lies on 
the cavity aperture and is the sum of the incident and reflected fields. Note that 
(4.37) contains the permeability of the exterior material {fire) while the corresponding 
function for the FE-BI formulation (3.25) contains the interior material (or cavity fill) 
permeability (p ri ). This is a consequence of applying Green’s first vector identity to 
either the exterior or interior volume regions, respectively. Thus, the same aperture 
excitation term may be used for the F E-ABC and FE-BI formulations. 

4.3.2 Material Term 

As previously mentioned, if a material is present in the exterior region, an ad- 
ditional volume and surface source terms need be included in (4.12). For example, 
these terms are necessary for a printed antenna with a radome covering. Such terms 
are a consequence of using a scattered field formulation in the exterior region. Typi- 
cally, a radome is a complex, inhomogeneous layered structure which exhibit desired 
transmission and reflection characteristics within a design bandwidth. Traditional 
method of moments formulations can be used to model such complex layered over- 
lays; however, such an approach requires an extremely complex and costly dyadic 
Green’s function or alternatively a full volume formulation, either of which is as- 
sociated with a fully populated matrix. Since the exterior region in this FE-ABC 
formulation is modeled using finite elements, radomes or any material covering may 
have arbitrary composition without the necessity of a complex Green’s function. In 
fact, the same algorithm may be used for both coated and uncoated bodies by simply 
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including the material source terms which are given next. 

The volume contribution may be recognized as the third term on the right-hand 
side of (4.11) 

f 7 = ~k 2 0 — -t T [ W, ■ E cyl dV (4.3S) 

[fir J JV‘ d 

This contribution is associated with the volumetric electric field and using the cylin- 
drical shell elements (3.9) presented previously, this term is given by 
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(4.39) 


where the subscripts indicate the orientation of the i th edge. It remains to determine 
the electric field anywhere in the exterior region which will require traditional modal 
analysis. 

An incident plane wave (also presented in (3.26)) is given by 
E { = 

= j/5 (sin 7 sin 4>< -f cos 7 cos 6{ cos </>,) + 4> (sin 7 cos </>,• — cos 7 cos 9{ sin <f>^ — 

z cos 7 sin 6i e i M' ,8in e< cos * t+z cos *] 

H' = Y 0 ^ p (sin 7 cos 0, cos <^, — cos 7 sin — 4> (sin 7 cos sin 0, + cos 7 cos 0,-^ — 
z sin 7 sin 6, e jko l psinB ' 

which is assumed to be incident on the cylinder from the direction and 7 

is the polarization angle where e = cos 7 # + sin 7 ^ is the electric field polarization. 
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In these expressions, the difference between the observation and incidence angles is 
denoted by <p, = 0 — o t . 

The electric field associated with the incident and reflected fields is determined 
by traditional modal analysis which enforces a vanishing tangential electric field on 
the surface of the cylinder (p = a). This field is given by 
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where the modal coefficients are given by 
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( 4 . 41 ) 


In these, J n (•) and Y n (•) are Bessel functions of the first and second kind, respec- 
tively. Primed expressions in (4.41) denote differentiation with respect to the argu- 
ment (e.g. J' n {0 = (C))- Note that the polarization angle ( 7 ) should not be 

confused with the modal coefficient ( 7 ^). For reference, the corresponding magnetic 
fields are given by 
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If (4.42) is evaluated on the surface of the cylinder (p = a). Q n = = 0 as required 

since the normal component of the magnetic field (or the tangential component of 
the electric field) must vanish on a metallic surface and (4.42) reduces to (3.28). 

These doubly infinite series may be converted into singly infinite series which are 
more amenable to numerical implementation 
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All of the formulae given above require numerical evaluation of the integrals which 
may be performed with a mid-point rule. At each integration point, the magnetic 
field is determined by truncating the infinite series in ( 4 . 43 ) so that the field value 
satisfies the convergence criterion 

yilK^n) - K»‘(n - 1)11 


H cyl — 


vW(n - Oil 


< 0.00001 


(4.45) 
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where u € {p.O.z}. Since these series exhibit poor convergence for large ratlins 
cylinders, the modal formulae (4.-13) are only suitable for relatively small radius 
cylinders (i.e. k 0 sin 6, a < 50). For larger radius cylinders, asymptotic expressions 
similar to (3.29) and (3.31) are more practical. 

When there is a discontinuity in permeability between material layers, a surface 
term must also be included on the right-had side of (4.11). That term is given by 
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and may be evaluated in a similar manner as the ABC surface term (4.37). After 
some manipulation, (4.46) is given by 
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assuming that the curved integration surface has h — p and that the cylindrical shell 
shape functions (3.9) are used for field representation. For an interface surface which 


has h — — the functional (4.46) is evaluated as 

rs d ^ w [ 1 1 1 r r ( z -J«) 


/f/ = jk 0 z 0 s ~^ — r r b{ ^-^H?{pAi,z)dpdz 

h Pre Pri J J z b ^ Pa P 

f? d = -jhzM- -1 [“ f n (p-p,)Hf{p,^,z)dpdz (4.48) 

t pre Pri \ z b ^ Pa 


while if h = <i>, (4.46) is given by 


iPb [ 1 1 1 f Zt [»» (z - Zi) 


iKz 0 -~f. — ~ — r r 

n Pre pri _ ** z b ^ Pa 


Hf (p,4 r ,z)dpdz 


/“ jKzM- -1 P f’\p-k) H fMr,z)ipdz (4.49) 

t Pre pri J ^ z b ^ Pa 

For an interface surface which has n = —z, the functional (4.46) is evaluated as 
ff/ = -jk.zjm f— - — 1 [*' r (4 _ *) (p, 4 , Zl ) dpd<t> 

P 0 [Pre Pri\ Ul J Pa V ' 

fii = jk 0 Z 0 y — — 1 f [ (p - pi)H c / {pA,z b ) pdpdcj) (4.50) 

0 t /ire Pri J -V 4po 
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while if ti = (4.46) is given by 


Sfi = -jKz/ ,Pt 


Q 


/a = 


5 t - 


--- rr(»- i.) p.o. : ,)d P do 
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The magnetic fields are computed using (4.44) -where the series are truncated accord- 
ing to (4.45). 

With the specification of the surface and volume source functionals, the FE-ABC 
methodology is now fully developed for a full range of typical applications. For 
radiation problems, since the total field is equal to the radiated field, the FE-ABC 
formulation is quite similar to the FE-BI formulation with the exception of replacing 
the exact mesh closure condition on the surface of the cylinder with an approximate 
condition at some distance above to the cylinder’s surface. However, for scattering 
analysis, in addition to the new mesh closure surface/condition, the excitation term 
may include radome contributions. Thus for coated antennas, scattering analysis 
requires additional volume and surface source terms which are a consequence of 
using a scattered field formulation in the exterior region. 



CHAPTER V 


Scattering 


One of the principal goals of this research was to conduct studies of the scat- 
tering by cavity-backed antennas recessed in a circular cylinder. Over the course of 
some forty years of radar cross section (RCS) investigations, the electromagnetics 
community has learned techniques which suppress traditionally dominant scattering 
contributions such as reflected and creeping waves associated with curved surfaces. 
Having reduced these bulk scattering features through the use of exotic materials 
and platform shaping, previously insignificant scattering contributions have become 
increasingly important. Scattering by an array of cavity-backed antenna elements is 
one such secondary contribution. 

In the past, most conformal antenna array designs were performed using either 
rigorous numerical formulations involving an infinite ground plane or approximate 
techniques such as the cavity model. Planar solutions are suitable for near broadside 
scattering provided the cylinder radius is large. However, since no circulating waves 
are supported by a metallic plane, bistatic scattering analysis in the penumbra and 
shadow regions of a cylinder cannot be obtained from planar results. Additionally, 
creeping wave coupling creates a significant scattered field on the rear-side of the 
cylinder even for backscatter measurements. Modern military aircraft demand low 
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observables over a full range of aspect angles necessitating the inclusion of creeping 
wave effects in any useful analysis technique. Approximate cavity models are use- 
ful for radiation analysis at resonance; however, they are unsuitable for scattering 
applications. 

The two hybrid finite element techniques presented in this thesis provide a means 
of performing, in the case of the FE-BI method, a rigorous analysis of scattering by 
cavity-backed antenna elements which are recessed in an infinite, metallic cylinder. 
The FE-ABC method provides an approximate but more flexible technique which 
may include the presence of a radome. In this chapter, the FE-BI method will be 
validated by comparison with data generated by different numerical models. This 
method will be shown to be highly accurate and, as implemented via a BiCG-FFT 
solver, highly efficient in terms of computational and memory demand. The FE- 
BI method will then be used to study the scattering behavior of typical conformal 
antenna elements. In particular, the effect of cavity size and curvature will be inves- 
tigated. The results will suggest useful design concepts which will aid RCS engineers 
in developing low observable antenna arrays. 

The FE-ABC formulation presented in this thesis will be compared to the re- 
sults generated by the FE-BI formulation. Their agreement will confirm, for the first 
time, that ABC termination schemes may be used above infinite, metallic structures. 
Previously, all FE-ABC applications involved closed, finite bodies. The FE-ABC 
approach offers considerable advantages over the FE-BI method since it may be effi- 
ciently used for computing the scattering by conformal antennas recessed in doubly 
curved platforms. Thus, it is essential that this technique be validated with a rig- 
orous solution, before its application to the necessarily more complex doubly curved 
structures. Additionally, the FE-ABC method provides a means of readily including 
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a complex, inhomogeneous radome which is inevitably placed above a conformal ele- 
ment. This approach also permits the modeling of a radar absorbing material (HAM) 
coating which is typically placed around a transparent radome. Such a RAM layer 
is used to suppress creeping waves and, therefore, it is an important feature in low 
observable designs and it is also useful in suppressing Electromagnetic Interference 
(EMI). 

5.1 FE-BI: Validation 

Having solved for the electric fields induced by an incident plane wave, the com- 
puted RCS data must be validated with known results. Available measured or com- 
puted data are rather scarce and as a consequence, we are forced to rely on limiting 
cases in order to validate this work. For example, as the radius of curvature de- 
creases, a cylindrical-rectangular cavity will approximate a planar-rectangular cav- 
ity. Another limiting case involves comparison of an elongated 3-D cavity with a 
corresponding 2-D cavity for normal incidence (0, = 90°). Finally, the infinite cylin- 
der results are compared with a finite Body of Revolution (BOR) model for oblique 
angles of incidence. 

The first validation effort for scattering by cavity-backed patch antennas relies 
on the fact that a small patch on a very large radius cylinder is quasi-planar and 
approximates rather well a planar-rectangular patch. For our test we chose as a 
reference a planar 3.678 cm x 2.75 cm patch residing on a 7.34 cm x 5.334 cm x 
0.1448 cm cavity filled with a dielectric having e r = 4. The equivalent patch on a 32.6 
cm cylinder is 6.46° x 2.75 cm residing on a 12.90° x 5.334 cm x 0.1448 cm cavity. 
At the operating frequency of 9.2 GHz, the cylinder has an electrical radius of 10A o . 
Figure 5.1 shows the results for the patch on a large radius cylinder with correspond- 
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Figure 5.1: Comparison of RCS for a planar patch (3.678 cm x 2.75 cm) residing on 
a 7.34 cm x 5.334 cm x 0.1448 cm cavity filled with t r — 4 dielectric and 
a corresponding quasi-planar patch on a large radius (10A o ) cylinder. 
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ing data for the planar cavity-backed patch. Clearly, the two RCS patterns are in 
excellent agreement, and although Figure 5.1 illustrates only monostatic scattering 
in the 0 = 0° plane, additional runs for normally incident monostatic scattering and 
various bistatic situations yield similar agreement. 

Comparisons may also be made for elongated cavities and 2-D moment method 
results. Long narrow cavities have very little axial interaction for principal plane 
{0 = 90°) incidence and therefore results based on this formulation should compare 
well with corresponding 2-D data. It is well known that the RCS of the 3-D scattering 
body of length L A 0 is related to the corresponding 2-D scattering of the same 
physical cross section via the relation 



Such a comparison is shown in Figure 5.2 for monostatic scattering by a 45° x 
5 A 0 x 0 .lA o cavity recessed in a cylinder with a radius of 1 A 0 for both principal 
polarizations. Once again the agreement between the two results is excellent, thus 
providing a partial validation of the formulation for highly curved geometries. Similar 
agreement has been observed for bistatic scattering in the 6 = 90° plane. 

The planar approximation eliminates the effects of curvature, which is of primary 
interest in this work, and the 2-D comparisons done above are only valid for normal 
incidence. To consider oblique incidence on a highly curved structure, comparisons 
with a Body of Revolution (BOR) code are made for wraparound cavities. Since 
the BOR code can only model finite structures, an infinite cylinder is simulated by 
coherently subtracting the far-zone fields of the finite structure without a cavity 
from similar data which includes the cavity. Such an procedure mimics common 
measurement practices and was found suitable for near normal incidence in the case 
of H-polarization (7 = 90°). This procedure is illustrated in Figure 5.3. An example 
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Figure 5.2: Comparison of 2-D moment method results and FE-BI results for a 45° x 
5A 0 x 0.1A o air-filled cavity which is recessed in a 1A 0 cylinder. 






74 


Far-Zone Subtraction Procedure: 



Ogive End-Cap °8 Jve End-Cap 

minus 




Figure 5.3: Far-zone subtraction scheme for the simulation of an infinite cylinder. 



calculation for the latter case is given in Figure 5.4 where a bistatic scattering pattern 
is presented in the d> = 0° plane due to a plane wave incident at (0, = 80°, o, = 0 1 ')- 
The BOR model was a 10A o long cylinder with 36° ogival end-caps. The wraparound 
groove had a height of 3A 0 and a depth of 0.1A„. Clearly, there is good agreement 
between the FE-BI results and data based on the BOR formulation. 

5.2 Scattering Studies 

The previous comparisons serve to validate the formulation. Having done so, it 
is instructive to examine the effect that curvature has on the scattering properties of 
cavity-backed patch antennas. Consider a 2 cm x 3 cm patch residing on a 5.0 cm x 
6.0 cm x 0.07874 cm cavity which is filled with a dielectric having t T = 2.17. Figures 
5.5 and 5.6 illustrate the behavior of this geometry as a function of frequency and 
curvature for E- and H-polarization, respectively. Evidently, the resonance behavior 
of this patch is sensitive to curvature for both principal polarizations. The frequency 
response for E-polarization is more sensitive to curvature since the radiating surface 
field component is parallel to the long side of the patch and cavity. If the patch and 
cavity were oriented so that the long side is parallel to the ^-direction, the response 
to H-polarization would exhibit greater sensitivity. Such an effect is important to low 
observable antenna designers since they want to operate the antenna in the region 
of lowest RCS. This low return region is a consequence of delicate cancellations due 
to the physical layout of the aperture. Such cancellations are not as complete for 
highly curved structures as they are for planar cavities. 

Conformal antenna designers often use wraparound antenna arrays to achieve om- 
nidirectional coverage. Two different configurations are typically used: a wraparound 
cavity where the cavity is filled with a single continuous collar of dielectric and dis- 
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Figure 5 . 4 : Comparison of the RCS computed via the FE-BI method and a BOR code 
for a 3A 0 x 0.1 \ D air-filled wraparound cavity excited by a H-polarized 
(7 = 90°) plane wave with oblique incidence ( &i = 80 = 0 °). 
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Figure 5.5: RCS frequency response for a 2 cm x 3 cm x patch residing in a 5 cm 
x 6 cm x 0.07874 cm cavity with t T = 2.17 as a function of curvature 
for E-polarization (7 = 0°). 
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Figure 5.6: RCS frequency response for a 2 cm x 3 cm x patch residing in a 5 cm 
x 6 cm x 0.07874 cm cavity with t r = 2.17 as a function of curvature 
for H-polarization (7 = 90°). 
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crete cavities symmetrically placed around the circumference of the cylinder. These 
two configurations are shown in Figure 5.7. Since near resonance, the radiation prop- 
erties of these two types of antennas are similar if all elements are uniformly driven, 
any RCS advantage which one might possess could govern the appropriate choice of 
arrays. Figure 5.8 compares the E-polarized monostatic scattering at 3 GHz in the 
0 = 90° plane for a continuous wraparound cavity and four discrete cavities; where 
the patches and cavities are identical to those used in the previous example. Not sur- 
prisingly, the continuous wraparound structure has a higher return due to coupling 
within the substrate. However, since in this case the scattered field is due to the 
^-component of the surface field ((^-directed magnetic currents), both cavities yield 
large scattered fields in the four directional lobes. Figure 5.9 is the corresponding 
comparison for H-polarization. In this case, the scattered field is attributed to the 
(^component of the surface fields (z-directed magnetic currents). Therefore, sub- 
strate modes diffract near the patch resulting in discrete lobes for the discrete array 
while creeping waves shed isotropically for the continuous wraparound cavity. Low 
observable designs will favor discrete cavity arrays over wraparound cavities since the 
scattering may be channeled in preferred directions and the overall scattering level 
is consistently lower. A final example is shown in Figure 5.10 where we observe that 
other than the expected higher scattering from the wraparound cavity, the scattering 
behavior of the two arrays in the 0 = 0° plane for H-polarization is very similar. 

5.3 FE-ABC: Validation 

In the previous section, the FE-BI formulation was used to investigate the role 
that both curvature and cavity dimensions play in determining the overall scattering 
properties of cavity-backed conformal antennas. The FE-BI method is extremely 
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(a) 



(b) 


Figure 5.7: Illustration of two types of wraparound arrays: (a) continuous cavity 
array; (b) discrete cavity array. 
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Figure 5.8: Comparison of E-polarized monostatic RCS at 3 GHz for a four patch 
array placed on a wraparound collar or in four discrete cavities. The 
patches and cavities are identical to the ones used in Figure 5.5. The 
observation plane is 0 = 90°. 
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Figure 5.9: Comparison of H-polarized monostatic RCS at 3 GHz for a four patch 
array placed on a wraparound collar or in four discrete cavities. The 
patches and cavities are identical to the ones used in Figure 5.6. The 
observation plane is 6 = 90°. 
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Figure 5.10: Comparison of H-polarized monostatic scattering at 3 GHz by a four 
patch array placed on a wraparound collar or in four discrete cavities. 
The patches and cavities are identical to the ones used in Figure 5.6. 
The observation plane is <f> = 0°. 
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powerful and efficient when implemented with a BiCG-FFT solver. However, as 
previously mentioned, it is not very flexible due to the uniform surface gridding re- 
quirement and limited number of available dyadic Green's functions. Rather, the 
FE-ABC method holds greater promise for impact since it can readily be used to 
model arbitrary patch and cavity shapes, coated and doubly curved conformal an- 
tennas. 

To date, the FE-ABC method has only been used to study the scattering by 
finite bodies. In this section, the FE-ABC method will be compared with rigorous 
FE-BI results for typical uncoated cavity-backed antennas which are recessed in an 
infinite metallic cylinder. Since the ABC surface must be displaced a distance from 
the cavity aperture in order to recover accurate results, the minimum displacement 
will be determined for computational domain minimization. Since the ABC absorbs 
outgoing waves, each principal polarization and typical incident angle directions will 
be examined. 

Consider a 2 cm x 3 cm patch antenna which is printed atop a metallic cavity 
encasing a 5 cm x 6 cm x 0.07874 cm substrate which has a dielectric constant of 
e T = 2.17 and scattering calculations are made at 3 GHz. The second-order ABC 
(4.25) is placed r\ 0 from the cavity aperture while the lateral walls (e.g. fi — ±<^ or 
h = ±£) of the ABC are placed 0.5A o from the cavity aperture. Experience suggests 
that 0.5A o lateral displacement is sufficiently large so that the only surface that need 
be considered is the normal ( h = p) surface which is placed at different distances, 
r, from the cavity aperture as shown in Figure 5.11. The bistatic scattering for 
normal incidence {(j) x — O°,0 t = 90°) and observation in the 0 — 90° plane due to an 
E-polarized plane wave is shown in Figure 5.12. To determine of the effectiveness 
of the ABC, Figure 5.13 illustrates a ten degree section of Figure 5.12 near the 
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backscatter region. Evidently, if the ABC surface is placed at least 0. 3 A 0 from the 
surface, the FE-ABC method recovers the FE-BI data to an acceptable accuracy. 
The corresponding patterns for H- polarization are shown in Figures 5.14 and 5.15, 
respectively. For this polarization, the FE-BI results are recovered when the ABC 
surface is only 0.1 A 0 from the cavity aperture. Such discrepancy in performance 
between the two principal polarizations illustrates the main drawback of the FE- 
ABC method. Since the ABC absorbs only the propagating waves, it must be placed 
beyond the region where near-fields are significant. However, the near-field/far- 
field boundary is problem and excitation dependent. Nevertheless, experience has 
shown that an ABC displacement of 0.3A o is a good rule-of-thumb for scattering 
computations. 

In the previous examples, the incident field was normal to the patch where the 
ABC is seen to work quite well. However, the performance of the ABC must be 
shown for grazing and rear-side incidence as well as normal. To do so, the same 
antenna element used previously, with the ABC surface displaced 0.3A o from the 
cavity aperture, is compared in Figures 5.16 and 5.17 with the FE-BI results for 
observation in the 6 = 90° plane. The FE-ABC method is seen to recover the 
FE-BI results even for incidence from the rear side of the cylinder. 


0 - 2 - - 



86 



Figure 5.11: Cavity-backed patch antenna with ABC mesh termination. 
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Figure 5.12: Convergence of the second-order ABC as a function of displacement 
from the cavity aperture for E-polarized bistatic scattering. The ref- 
erence data is provided by a rigorous FE-BI formulation for the same 
cavity-backed antenna. 







89 


CQ 

T3 




GO 

u 

06 

u 

• y< 

cS 

4-> 

cn 

OQ 



Angle (<|>) [deg] 


Figure 5.14: Convergence of the second-order ABC as a function of displacement 
from the cavity aperture for H-polarized bistatic scattering. The ref- 
erence data is provided by a rigorous FE-BI formulation for the same 
cavity-backed antenna. 
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Figure 5.15: Magnified view of Figure 5.14. 
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Figure 5.16: Monostatic scattering by a 2 cm x 3 cm x patch residing in a 5 cm x 
6 cm x 0.07874 cm cavity with t r = 2.17 for E-polarization (7 = 0°). 







CHAPTER VI 


Radiation 


In the previous chapter, the scattering behavior of cavity-backed antennas which 
are recessed in a metallic cylinder were investigated. Both the FE-BI and FE-ABC 
methods were used and found to be accurate as compared to available reference data. 
In this chapter, the radiation parameters of conformal patch antennas mounted on 
singly curved surfaces are examined. In particular, the antenna pattern and input 
impedance of such patch antennas are investigated. 

Approximate analysis methods such as the cavity model have been used to design 
conformal antennas. Initially, planar solutions [3] were obtained which have been 
found to match measured antenna patterns for operation at resonant frequencies. 
Since the cavity model assumes single mode fields below the patch, such agreement is 
understandable. However, in contrast, input impedance calculations using the cavity 
model are not accurate with respect to either input resistance or resonant frequency 
localization. The reason for this inaccuracy is the fact that a typical probe feed 
generally excites multiple modes and the fact that, in reality, the field is not confined 
to the region immediately below the patch. This latter bias is partially compensated 
by utilizing effective permittivities rather than measured permittivities. Typically, 
with this empirical modification, the resonant frequency localization is acceptable; 
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however, the input resistance calculations remain inaccurate in comparison with 
measured data. 

Although the cavity model is limited to resonant patches and its input impedance 
predictions are inaccurate, the cavity model has been successfully extended to rectangular- 
cylindrical patch antennas [4, 5, 53]. Such extensions to the cavity model have been 
used to explore the effect of curvature on both the radiation pattern and the input 
impedance of curved patches. Once again, the calculated patterns compare favor- 
ably with measurements [5, 54] while the input impedance calculations are rather 
poor. Use of an effective permittivity improves the resonant frequency localization; 
however, the input resistance is still unacceptably inaccurate for design use. 

Patch antennas fed by a probe have a very small bandwidth. In an effort to 
analyze printed broadband antennas such as spirals, rigorous numerical techniques 
based on integral equation formulations have been presented in the literature [6]- [9]. 
Since the metallization pattern is assumed to be printed atop a dielectric coated 
groundplane, all of these formulations employ a complex dyadic Green’s function. 
Typically, a rather simplistic feed model for a current probe is utilized and although 
the radiation pattern of such a model agrees with measurement, the input impedance 
once again has unacceptable error. Recently, Aberle, Pozar and Birtcher [55] have 
investigated improved feed models which include a finite probe diameter and an 
impressed voltage source rather than a current source. Although most numerical 
formulations require a planar structure, these integral equation techniques have also 
been extended to patches printed on a coated cylinder [11, 12]. Even though such 
formulations have proven useful for antenna analysis and design, the assumption of 
a coated, infinite substrate is not realistic. In practice, spurious coupling between 
elements within an array will occur on such a structure and as a consequence, each 
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antenna element is typically placed within a substrate filled metallic cavity. Aberle 
recently developed an integral equation formulation for planar cavity-backed an- 
tenna elements [13]. However, this formulation suffers, as do most integral equation 
techniques, from a very large computational and memory burden due to the fully 
populated matrix and complex dyadic Green's function. If uniform zoning is used 
on the patch surface, the resulting matrix will be block Toeplitz and accordingly it 
may be efficiently solved using the BiCG-FFT solver employed in this work. How- 
ever, such an implementation would still require a highly complex dyadic Green’s 
function. In addition, inhomogeneous substrates or radomes may not be efficiently 
solved using a surface integral equation. 

In response to the inefficiencies and limitations of the integral equation formu- 
lations, the FE-BI method was introduced by Jin and Volakis [17]. This approach 
has proven quite successful for planar antenna radiation pattern calculations. As 
discussed earlier, the inaccurate input impedance data computed by that code may 
be remedied by using an accurate feed model which is currently under investigation 
at the Radiation Laboratory. In this thesis, the FE-BI formulation is extended to 
circular cylinder platforms. 

In this chapter, the two finite element formulations, the rigorous FE-BI method 
and the approximate FE-ABC method, will be utilized to study cylindrical- rectan- 
gular patch antenna radiation. First, the FE-BI method will be validated for antenna 
pattern applications by comparison with measured data. The input impedance cal- 
culations will be seen to suffer from the same inaccuracies as the planar implemen- 
tation due to the use of a similar simplistic probe. After validation, we will examine 
the effect of curvature on the pattern and the input impedance of the cylindrical- 
rectangular patches. An important class of antennas which wrap around the cir- 
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cumference of the platform cannot be studied using planar models. Such antennas 
are used for communications and guidance control since omnidirectional coverage is 
essential. The proposed FE-BI has the capability of examining three typical types of 
wraparound antennas: discrete cavity, continuous cavity and single collar antennas. 
Since all three types of antennas may be analyzed by the FE-BI method, comparisons 
may be made between the three in terms of coverage and mutual coupling between 
elements. Finally, as was done in Chapter V for scattering calculations, the FE-ABC 
method will be compared to the FE-BI method and a minimal ABC displacement 
will be established. 

6.1 FE-BI: Validation 

Two types of patch antenna elements are investigated and these are shown in 
Figure 6.1. Each patch is aa° x b in size where a denotes the radius of the cylinder 
and a is the subtended angle of the patch. A patch whose radiating side walls 
form constant ^-surfaces is termed an axially polarized patch and is fed at <f>j = 
Circumferentially (or azimuthally) polarized patches have radiating walls forming 
constant z-surfaces and are typically fed at zj = Observation in the 6 = 90° plane 
is the E-plane for circumferentially polarized patches and the H-plane for axially 
polarized elements. The terminology originates with the cavity model for patch 
antennas. We will now characterize a typical cavity-backed patch antenna. 

Several computed and measured antenna patterns have been published for patches 
printed on a coated cylinder. One such patch, which is 3.5 cm x 3.5 cm, was used 
by Sohtell [54] to compare the accuracy of the cavity model [5] to a surface current 
integral equation [12]. The measured data was taken at 2.615 GHz for a metallic 
cylinder which was 63.5 cm long and had a radius of 14.95 cm. The cylinder was 
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Figure 6.1: Illustration of (a) a circumferentially polarized patch element; and (b) 
an axially polarized patch element. The radius of the cylinder is denoted 
by a. 
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coated with a 0.3175 cm uniform dielectric having relative permittivity of ( r = 2.3.2. 
Data was taken for —180° < o < 180° in the 0 = 90° plane corresponding to the 11- 
plane for circumferentially polarized elements and the H-plane for axially polarized 
patches. Figure 6.2 compares these measured patterns with data generated using 
the FE-BI formulation for an identical patch placed within a 360° x 7 cm cavity. 
This wraparound cavity best simulated the measured coated cavity. Note that the 
H-plane patterns are symmetric due to the symmetric placement of the feed, whereas 
the E-plane patterns are not symmetric. The placement of the feed was not specified 
in [54]; however, the agreement for the E-plane pattern shown in Figure 6.2 indicates 
that the position used in the FE-BI model (a</>^ = -1 cm) is reasonable. The feed 
was placed at zj = -1 cm for the axially polarized (H-plane) case. 

In addition to single patches, the FE-BI formulation may be used to design mi- 
crostrip arrays. Such an approach includes mutual coupling between elements which 
is ignored by the cavity model. Furthermore, the FE-BI formulation consumes less 
computational resources than a comparable integral equation formulation due to the 
sparsity of the FE matrix. The H-plane pattern of a four element array was measured 
to gauge the accuracy of the FE-BI approach. Each element is 2 cm x 3 cm and 
placed within a 5 cm x 6 cm x 0.07874 cm cavity which is filled with a dielectric 
having t r = 2.17. The cylinder is 91.44 cm long and has a radius of 15.24 cm. The 
cavities are placed symmetrically around the cylinder (e.g. a patch is centered at 0°, 
90°, 180° and 270°). Only the patch centered at 0° was excited while the remaining 
patches were terminated with a 5014 load. The driving patch is axially polarized and 
the feed is located at zj = —0.375 cm. Figure 6.3 illustrates the excellent agreement 
between the FE-BI formulation and the measured H-plane data. 
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Figure 6.2: Comparison of measured [54] and computed data for a circumferentially 
polarized element (E-plane) and an axially polarized element (H-plane). 
The 3.5 cm x 3.5 cm antenna was printed on a cylinder with a radius 
of a = 14.95 cm and a 0.3175 cm coating having e r = 2.32. The probe 
feed was place at (a 4>/,Zf) = (-1.0, 0.0) for the circumferentially polarized 
patch and at (a <f>j,zj) = (0.0, -1.0) for the axially polarized antenna. 
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Figure 6.3: H-plane pattern for a four element patch array. Each of the patches is 2 
cm x 3 cm and are placed symmetrically around the cylinder. Only the 
patch centered at 0° is fed while the other patches are terminated with 
500 loads. 
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6.2 Radiation Studies 

In Chapter V, discrete cavity arrays were found to have a significantly lower 
radar cross section compared to a continuous wraparound array. Thus, the size of 
the cavity had a significant effect on the scattering properties of the array. The 
axially and circumferentially polarized antennas used by Sohtell [54] were placed 
within cavities which were 7 cm high and approximately 30°, 50°, 90°, 180°, 270° 
or 360° in angular extent. Figure 6.4 illustrates that azimuthal cavity size has little 
effect on the radiation pattern for a circumferentially polarized element. A similar 
comparison for the axially polarized patch is shown in Figure 6.5. The back lobe 
of the antenna (near <j> = 180°) is very small for cavities less than 180° in extent 
but increases for larger cavities. For cavities which lie on the forward face of the 
cylinder, the substrate modes diffract off the cavity walls; an effect which has little 
influence on the main lobe of the pattern. However, for wraparound cavities and 
cavities which extend into the back side of the cylinder, the substrate modes shed 
like creeping waves giving rise to the back lobe. 

Having established the effect of cavity size on the antenna patterns, it is instruc- 
tive to gauge the effect that curvature has on the resonance behavior (or gain) of 
patch antennas. The two antennas were placed in 14 cm x 14 cm cavities recessed 
in cylinders with increasing radius. The frequency was allowed to vary from 2.4 GHz 
to 2.7 GHz and the peak radiated power (|i? 0 ( ^| 2 ) was recorded at each frequency. 
Figure 6.6 illustrates that the resonance frequency increases with increasing curva- 
ture for a circumferentially polarized antenna, but the maximum radiated power 
is similar regardless of element curvature. Note in the cavity model, the radiating 
edges for a circumferentially polarized patch are the azimuthal walls of the cavity 
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Figure 6.4: Effect of cavity size on the E-plane radiation pattern of a circumferen 
tially polarized patch antenna. 
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Figure 6.5: Effect of cavity size on the H-plane radiation pattern of an axially polar- 
ized patch antenna. 
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Figure 6.6: Resonance behavior of a circumferentially polarized patch antenna for 
various cylinder radii. 
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(see Figure 6.1) which have a constant separation regardless of the cylinder radius. 
However, the axially polarized patch has decreasing resonant gain with increasing 
curvature as shown in Figure 6.7. For this patch, radiation is attributed to the 
axial magnetic walls of the cavity model which have increasing angular separation 
with decreasing curvature. These walls radiate strongly away from the pattern peak 
(d> = 0°). Accordingly, the gain of an axially polarized antenna decreases with in- 
creasing curvature. The radiation pattern of a circumferentially polarized antenna 
is largely unaffected by curvature as shown in Figure 6.8 when excited at a resonant 
frequency. However, the radiation pattern of the axially polarized antenna broad- 
ens as the curvature increases which is illustrated in Figure 6.9. Once again, both 
relationships are readily explained by considering the effect curvature has on the 
orientation of the cavity model radiating walls. 

In addition to the gain and pattern of an antenna, designers require the input 
impedance for matching purposes. For the antenna examined above (in a 14 cm x 14 
cm cavity), the input impedance was calculated from 2.4 GHz to 2.7 GHz for various 
cylinder radii. Figure 6.10 illustrates that the input impedance of a circumferentially 
polarized patch antenna is not affected by curvature while Figure 6.11 shows that 
increased curvature reduces the input impedance of an axially polarized patch. This 
observation agrees the the results reported by Luk et. al. [53]. 

An important application of conformal antennas on drones is for guidance control. 
An omnidirectional pattern is required for this application since the presence of a null 
would have a disastrous effect on the drone’s guidance capabilities. Omnidirectional 
coverage may be obtained by placing a wraparound array on the cylinder such as the 
ones shown in Figure 5.7. Consider the axially fed 2 cm x 3 cm patch used previ- 
ously. Figure 6.12 illustrates the coverage obtained by various discrete wraparound 



I 


QQ 

T3 


O 

ON 

II 

CD 


O 

II 

■e- 


<u 

£ 

o 

CLh 

• 4 — • 

.2 

otf 



Figure 6.7: Resonance behavior of an axially p 
cylinder radii. 



(f) [GHz] 


patch antenna for various 




Radiated Power [dB 


107 



Angle (4>) [deg] 


Figure 6.8: Variation of the radiation pattern shape with respect to curvature for a 
circumferentially polarized antenna. 
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Figure 6.10: Input impedance of a circumferentially polarized patch antenna for var- 
ious cylinder radii. The frequency range was 2.4 GHz to 2.7 GHz and 
the cavity size was 14 cm x 14 cm. 
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Figure 6.12: Coverage of various uniformly fed discrete wraparound arrays. Each 
patch is 2 cm x 3 cm and the feed point is at (<f>j = <f> C enten z f = —0.375 
cm) where </> center denotes the azimuthal center of each patch. The array 
is operated at 3 GHz. 





arrays (see Figure 5.7a) operated at 3 GHz which corresponds to a cylinder radius of 
1.52789A 0 . Evidently, as elements are added to the array, phase interference causes 
side lobes to develop which are seen as ripples in the pattern. The geodesic distance 
between each element center is: 4.8A 0 , 2.4A 0 , 1.2A 0 and 0.6A o for two. four, eight 
and sixteen patches, respectively. When sixteen elements are used, the inter-clement 
spacing is sufficiently small as to suppress side lobes and the desired omnidirectional 
pattern is obtained. Generally, the required number of radiators will scale with the 
radius of the cylinder. Figure 6.13 shows similar results for a continuous wraparound 
array (see Figure 5.7b). 

Another common type of wraparound antenna utilizes a single, continuous metal- 
lic collar as a radiating element. Figure 6.14 illustrates the pattern for such an an- 
tenna when uniformly excited probe feeds are place symmetrically around the collar 
with Zf — —0.375 cm. The inter-feed spacing is identical to the center element spac- 
ing used in the previous examples. Evidently, as the number of feeds increases, the 
currents supported by the collar become more uniform and for the case of sixteen 
feeds, the pattern is omnidirectional. However, notice that for this case, even when 
only eight feeds are used, ail dB ripple is obtained. 

6.3 FE-ABC: Validation 

As previously mentioned in Chapter V, the FE-ABC method is considerably 
more flexible than the FE-BI approach. However, the FE-ABC formulation requires 
a computational domain which extends some distance beyond the cavity aperture. 
In Chapter V, a ABC surface displacement of 0.3A o was shown to be sufficient for 
most calculations. In this section, similar comparisons are made for both radiation 
pattern and input impedance calculations. 



Normalized Gain [dB] 


113 



- One patch 
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- Four patches 
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Sixteen patches 


Figure 6.13: Coverage of various uniformly fed continuous wraparound arrays. Each 
patch is 2 cm x 3 cm and the feed point is at (<f>j = 4> C enteT,Zf — —0.375 
cm) where 0 cen * cr denotes the azimuthal center of each patch. The array 
is operated at 3 GHz. 
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Figure 6.14: Coverage of various uniformly fed collar wraparound antenna. The 
metallic collar has a height of 3 cm and is fed at symmetric points 
around the cylinder with zj = —0.375 cm. The antenna is operated at 
3 < 
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Consider a 2 cm x 3 cm patch antenna which was printed atop a metallic cavity 
encasing a 5 cm x 6 cm x 0.07874 cm substrate having a dielectric constant of t T 
= 2.17. This is the same antenna used for the FE-ABC validation in Chapter V 
for scattering calculations. The second order ABC was place t \ 0 from the cavity 
aperture while the lateral walls (e.g. n = ±0 or n = ±£) of the ABC were placed 
0.5A o from the cavity aperture. The H-plane antenna pattern for an axially polarized 
patch is shown in Figure 6.15 where the probe feed is placed at 4>j = 0° and zj = 
—0.375 cm and the operating frequency is 3 GHz. Figure 6.16 shows the main 
lobe peak value of this antenna for various ABC surface displacements. In this 
example, the radiation pattern calculated using the FE-ABC formulation was seen 
to be in good agreement with the pattern computed using the more rigorous FE- 
BI approach. Similar comparisons for a circumferentially polarized antenna element 
(e.g. 3 cm x 2 cm patch within a 6 cm x 5 cm x 0.07874 cm cavity with a radial 
feed at 4>j = — 1.41°,z/ = 0 cm) are given in Figures 6.17 and 6.18. Evidently, 
placement of the ABC 0.3A o away from the cylinder’s surface is sufficiently accurate 
for patterns obtained from circumferentially polarized as well as an axially polarized 
patch antennas. However, as illustrated in Figure 6.19, the FE-ABC method is 
not sufficiently accurate for input impedance calculations. For the antenna used 
in this comparison, the FE-ABC solution produces an erroneous value even when 
the ABC surface is 0.4A o from the aperture. Evidently, the FE-ABC method does 
not accurately model the fields below the patch radiator unless the ABC surface is 
placed prohibitively far from the aperture. Even a very large ABC displacement 
will be inaccurate since the increasing computational error associated with a very 
large mesh will dominate the decreasing approximation error on the ABC surface. 
However, as shown in the previous example, far-field calculations are more forgiving 
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Figure 6.15: Convergence of the second-order ABC as a function of displacement 
from the cavity aperture for axially polarized patch. The reference data 
is provided by a rigorous FE-BI formulation for the same cavity-backed 
antenna. 




Radiation Pattern fdB 





Radiation Pattern [dB] 


118 



Angle (<(>) [deg] 


Figure 6.17: Convergence of the second-order ABC as a function of displacement 
from the cavity aperture for circumferentially polarized patch. The 
reference data is provided by a rigorous FE-BI formulation for the same 
cavity-backed antenna. 
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Figure 6.19: Input impedance calculations as a function of ABC surface displacement 
from the cavity aperture. The frequency range was 3.0 to 3.2 GHz with 
data taken every 5 MHz. 
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and hence the FE-ABC method is useful for radiation pattern calculations. 



CHAPTER VII 


Conclusion 


7.1 Summary 

In this thesis, a new approach for analyzing cavity-backed cylindrically confor- 
mal antenna arrays was presented. Previously, such antennas were designed using ap- 
proximate planar cavity models or rigorous integral equation formulations. Although 
both of these techniques have been extended to coated cylinders, planar solutions 
were typically used even for curved platform designs due to their simplicity. Both 
techniques work reasonably well for resonant antenna pattern calculations. However, 
for scattering or input impedance calculations, the integral equation method has high 
computational demands whereas the cavity model is too approximate. Additionally, 
both methods typically involve metallic patches printed atop a coated groundplane. 
A cavity-backed integral equation formulation is available; however, the calculation 
of the internal dyadic Green’s function limits its utility. 

Due to our past favorable experience with the application of FE-BI formulations 
to planar cavity-backed planar arrays, in this work the FE-BI method was applied to 
curved conformal arrays. A critical feature of the formulation is the use of the itera- 
tive BiCG-FFT solver which greatly reduces the memory and computational burden 
as compared to a direct implementation. Thus by using uniform surface gridding, 
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impressive computation speeds have been obtained at the expense of grid flexibility. 
Another important aspect of this work is the use of creeping wave expressions for the 
calculation of the cylindrical dyadic Green's function. The traditional modal series 
expressions of this Green’s function imposes a prohibitive computational burden on 
the boundary integral calculation while the creeping wave series expressions are both 
efficient and accurate when used for a large radius cylinder. 

The FE-BI formulation has been successfully used to study the scattering be- 
havior of conformal antennas. Discrete cavity arrays have been shown to have a 
lower overall RCS as compared to equivalent continuous wraparound arrays. The 
primary reason is that the cavities suppress creeping and substrate wave coupling. 
Such waves diffract into preferred directional lobes which allow the designers to chan- 
nel the scattered field away from the source. Conversely, wraparound collar arrays 
radiate isotropically resulting in a more uniform pattern. Therefore, low observable 
antenna designers would prefer discrete arrays over continuous arrays due to the 
increased control available with the former conformal elements. 

The effect of curvature on scattering resonance was also examined. Evidently, 
for either principal polarization, the resonant scattering is not affected by curvature. 
However, depending on the orientation of the patch, the null which succeeds the 
resonant peak is curvature dependent. Highly curved patches have a more shallow 
null than planar or nearly planar patches. RCS engineers would like to place such 
a null in an expected threat frequency and accordingly, they would prefer a planar 
patches in such circumstances. 

In addition to scattering studies, the FE-BI formulation was used for radiation 
analysis. Although this implementation is sufficient for radiation pattern calcu- 
lations, as was the case with the planar formulation, the infinitesimally thin probe 
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feed was found to be too simplistic for accurate input impedance computations. Nev- 
ertheless, since this feed is sufficient for resonant frequency localization, the FE-BI 
method was used to study the effect of curvature on the resonant frequency of curved 
patch antennas. Evidently, the resonant frequency of an axially polarized element 
is sensitive to curvature while a circumferentially polarized element is insensitive to 
curvature. Likewise, the co-polarized pattern associated with the axially polarized 
antenna exhibited a marked increase in the back lobe while the pattern due to the 
corresponding circumferentially polarized element was essentially unchanged by cur- 
vature. Also, the axially polarized element’s co-polarized pattern is dramatically 
effected by cavity size while the corresponding circumferentially polarized pattern is 
barely altered by increasing curvature. Likewise, the input impedance of an axially 
polarized element is dependent on the curvature of the substrate while the circum- 
ferentially polarized element has a input impedance which is insensitive to curvature. 

Another important aspect of this thesis was the application of the FE-ABC 
method for radiation and scattering by cavity-backed antennas on an infinite struc- 
ture. In this, a mixed total/scattered field formulation was introduced which substan- 
tially reduced the likelihood of error propagation and decreased the computational 
effort of the implementation. New second-order conformal ABCs [25] were used in 
this FE-ABC method which greatly reduced the computational domain. An ABC 
displacement of only 0.3 A is usually sufficient for both scattering and radiation pat- 
tern calculations. However, no reasonable displacement was found which yielded 
accurate input impedance data. The FE-ABC method was found to be accurate for 
both bistatic and backscatter calculations even for incidence in the rear side of the 
cylinder. Thus, the presence of the ABC does not appear to alter the creeping wave 
interactions on the surface of the cylinder. Both principal polarizations were used 



for the scattering study as well as axial and circumferentially polarized radiating an- 
tenna elements. We expect that the FE-ABC method will be accurate for arbitrary 
incidence and polarization. 

7.2 Future Work 

Within the framework of this thesis, two obvious avenues of research are sug- 
gested: development of an improved feed model and application of the FE-ABC 
method to doubly curved and coated conformal antennas. Historically, the develop- 
ment of feed models for finite element applications has paralleled the similar develop- 
ment of feeds for integral equation formulations. Both methods initially employed a 
simple infinitesimally thin wire probe. Integral equations then used increasingly more 
advanced feeds until the current state-of-the-art, the base voltage gap. Such a voltage 
gap source was used by Aberle [55] where the current supported by the metal wire is 
determined via an integral equation. The finite thickness of the wire is simulated by 
placing the current expansion functions at the center of the wire while testing at the 
outer radius. Unfortunately, such an approach is not readily implemented in a finite 
element formulation though advanced FE feed models are desirable. For example, 
Gong and Volakis meshed a coaxial cable feed using tetrahedral elements to a point 
somewhat below the cavity where the lowest order coaxial mode was used as a coax- 
ial aperture source. This technique, although rigorous, yielded poorly conditioned 
matrices and imposed a prohibitively expensive computational and memory burden. 
A different approach is to apply the lowest order coaxial mode across the aperture of 
the cable at the base of the cavity itself. Such a model will work quite well at the res- 
onant frequency since in a matched condition, the only coaxial mode present at this 
aperture will be of the lowest order. However, for operation away from resonance, 
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the mismatch at the cable aperture would require higher-order coaxial modes. A key 
aspect of this feed model is to use a voltage continuity condition to couple the cable 
aperture fields to the interior fields. Thus, the challenge is in matching the voltage 
due to higher order modes with the potential within the cavity. If however a field 
continuity condition is used, mixed elements must be combined perhaps through the 
use of a boundary integral. 

Although this work is useful to the design community in and of itself, its main 
purpose was to provide a key step in the evolution of conformal antenna simulation 
from planar structures to doubly curved and coated platforms. The FE-BI technique 
was used to validate the FE-ABC method for singly curved antennas. Additionally, 
the FE-ABC method was seen to be accurate for scattering and radiation pattern 
calculations by cavity-backed structures which are recessed in an infinite metallic 
surface. With confidence, similar formulations may be developed for doubly curved 
and coated antenna elements. To do so, new distorted brick elements should be de- 
veloped. Such elements will not be divergence free and hence the method will require 
a penalty term; however, use of distorted bricks will allow custom mesh generators 
to be used which would alleviate the need to utilize a sophisticated (and expensive) 
general purpose mesh generator such as SDRC-IDEAS. The FE-ABC method should 
be used rather than the FE-BI approach since a doubly curved surface does not al- 
low use of a BiCG-FFT solver and hence the computational and memory burden 
of a FE-BI method is prohibitive. Additionally, the FE-ABC approach allows the 
introduction of complex, inhomogeneous material radomes. For such coated geome- 
tries, the equivalence principle will require two vector currents on the computational 
domain border for far-zone field computations. The reflected fields will need be 
determined for a doubly curved surface via high frequency techniques and efficient 
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on-surface and far-zone formulae will have to be developed. 
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APPENDIX A 
Fock Functions 


The asymptotic form of the dyadic Green’s function with observation both on the 
surface of the cylinder and in the far field involves Fock functions. These have been 
extensively studied and tabulated by Logan [40]. In this appendix, the numerical 
evaluation of these functions is performed for both small arguments and large argu- 
ments. 

The on-surface Fock functions used in this work are 

* ■ Wif 
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where w 2 (t) denotes Airy functions of the second kind and w' 2 (t ) is its derivative 
with respect to r. For small arguments (£ < 0.6), the asymptotic expansion of (A.l) 
is given by 
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while a rapidly converging residue series is used for £ > 0.6 
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Zeros of the u^(t) and ir^r) 

i 

T n 

= i TV, j £ J 3 

and T n = |r^|f 

n 

l r »l 

l r nl 

i 

2.33811 

1.011879 

2 

4.08795 

3.24819 

3 

5.52056 

4.82010 

4 

6.78661 

6.16331 

5 

7.94413 

7.37218 

6 

9.02265 

8.48849 

7 

10.0402 

9.53545 

8 

11.0085 

10.5277 

9 

11.9300 

11.4751 

10 

12.8288 

12.3848 


Table A.l: Zeros of the Airy Function 

10 

u(f) ~ =2 e ^^ l E( r n)" 1 e“^ Tn (A. 3) 

n= 1 

where r n and are zeros of w 2 (r) and u4(t), respectively. Those zeros are given in 
the Table A.l 

The far-zone Fock functions are given by 
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where u’i(r) denotes Airy functions of the first kind. u’,(t) is its derivative with 
respect to r and the appropriate integration contour (T) is given by Logan [40 j. 1 wo 
of the functions required in this work, g (0) (£) and <7 (1) (£), ma . v be calculated using 
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with constant k = e _ - ? W 6 and the coefficients for (A.5) and (A.6) are given in Table 
A. 2. In a similar manner, /(°)(£) may be computed using 
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which utilizes the constants in Table A. 3 
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Constants for (A. 5) and (A. 6) 

m 

c(m) 

a' (m) 

At(m) 

i 

0.7473831 

1.01879297 

0.5356566 

2 

-0.6862081 

3.2481975 

-0.41901548 

3 

-2.9495325 

4.82009921 

0.38040647 

4 

-3.4827075 

6.16330736 

-0.35790794 

5 

8.9378967 

7.37217726 

0.34230124 

6 

56.1946214 

8.48848673 

-0.33047623 

7 


9.53544905 

0.32102229 

8 


10.52766040 

-0.31318539 

9 


11.47505663 

0.30651729 

10 


12.38478837 

-0.30073083 


Table A. 2: Constants for <^°)(£) and g^(£). 
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Constants for (A. 7) 

m 

C (m) 

a(ra) 

Ai^m) 

i 

1.146730417 

2.33810741 

0.70121082 

2 

0.86284558 

4.08794944 

-0.80311137 

3 

-2.0192636 

5.52055983 

0.86520403 

4 

-9.977776 

6.78670809 

-0.91085074 

5 

-14.59904 

7.94413359 

0.94733571 

6 

49.0751 

9.02265085 

-0.97792281 

7 


10.04017434 

1.00437012 

8 


11.00852430 

-1.02773869 

9 


11.93601556 

1.04872065 

10 


12.82877675 

-1.06779386 


Table A. 3: Constants for /^(0* 
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APPENDIX B 


Iterative Solvers 


Iterative solvers are useful for large sparse linear systems since they typically require 
considerably less memory as compared to direct solvers such as LU decomposition. 
Two iterative solvers are used in this work: the biconjugate gradient and conjugate 
gradient squared algorithms. The former is used for the symmetric systems which 
are a result of using the FE-BI method or the FE-ABC method in conjunction with 
the symmetric ABC (4.32). The latter is used for the FE-ABC formulation which 
employs an asymmetric ABC such as (4.25) or (4.27) without any artificial sym- 
metry. This appendix will present both algorithms beginning with the biconjugate 
gradient. Note that the algorithms presented are do not employ any preconditioning. 
Preconditioned versions are readily found in the literature (see for example [52]). 

B.l BiCG Algorithm 

The biconjugate gradient (BiCG) algorithm is one member of a family of iterative 
solvers which have proven useful in computational electromagnetics [45]. The BiCG 
unlike the conjugate gradient (CG) method has the advantage of utilizing only one 
matrix-vector product in its symmetric implementation. Although the convergence 
characteristics of the BiCG algorithm are erratic (see for example [45]), it often 
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converges in fewer iterations than the CG algorithm. This section presents the BiCG 
algorithm appropriate for use with symmetric matrices [46]. 

Consider the system 

Ml W = W (B.l) 

where [/l] is a symmetric matrix, {x} is the unknown data vector and {6} is the 
excitation data vector. The BiCG pseudocode follows (assuming an initial guess 
{x}j which may be {0}): 

Initialize: 

Mi = M-MKK 
Mi = W-MH 


Iterate: 



<m„,m:> 
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< [A] M n , M^ > 


Mn+1 = 

{*}„ + «» M„ 


Mn+1 = 

M„ - fl » l A ] Mn 




^ { r }n+l >{ r )n+l ^ 


— 

< M„>Mn > 


m = 

M„+i + c„ {p} n 

(B.2) 

where the Euclidean norm is given by 



M 

< {x} , {6} > = x [ m ]^*[ m ] 

m— 1 

(B.3) 

In (B.2) the subscripts refer to the iteration and the asterisk denotes complex con- 


jugation. It is necessary to terminate the iteration when the error present in the 
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solution is suitably small. Although the error cannot be directly measured (other- 
wise there would be no need to solve the system), the error does control the correct ion 
factor between each iteration. When the error is large, the correction applied to the 
new guess is large while a small error results in only a small change in the result from 
one iteration to the next. Many termination criteria have been used in the past. One 
of the most popular is 


< Wn+llMn±l > < 

<{b},{b}> 


(B.4) 


where e < 0.01 is a typical acceptable tolerance. As can be expected, the tolerance 
may need to be tightened or relaxed depending on the problem at hand and the de- 
sired accuracy of the result. Note that t should be kept small for antenna impedance 
computations but can be relaxed for scattering and radiation pattern calculations. 


B.2 CGS Algorithm 


The conjugate gradient squared (CGS) solver, like the BiCG, does not guarantee 
convergence. Unlike the symmetric version of the BiCG algorithm, the CGS requires 
two vector-matrix multiplies; however, the CSG solver can be used for either sym- 
metric or asymmetric matrices. As is the case with the BiCG solver, the convergence 
of the CGS algorithm is erratic; however, it often requires fewer iterations than the 
either the BiCG or the CG algorithm. This section presents the CGS algorithm 
appropriate for use with either symmetric or asymmetric matrices [52]. 

Consider the system 

[A]{*} = {b} (B.5) 

where [A] is a symmetric matrix, {x} is the unknown data vector and {6} is the 
excitation data vector. The CGS pseudocode follows (assuming an initial guess {x}j 
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which may be {0}): 


Initialize: 


Mi 

Mi 

H 

i 

-O ^ 

II II 

(B.6) 

Mi 

= {1} 

(B.7) 

Mi 

= {1} 

(B.8) 

p i 

= 0 

(B-9) 


(B.10) 


Iterate: 

Pn = < Ml , Mn-1 > 

/*=-£=- 
Pn- 1 

M = Mn-1 +P M„-l 
M« = M +/?({<?}„_! + /?M„-l) 

W = MM„ 

“ = < mTTm > 

M„ = {w}-Or{u} 

M = M + M n 

{*}„ = {x} n _!+a{u;} 

M„ = Mn-1 “ a MM} (B.ll) 


The same termination criterion as was used for the BiCG algorithm is used for the 
CGS algorithm. 
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APPENDIX C 


FFT-based Matrix- Vector Product 


The numerical solution of integral equations (IE) (or boundary integrals) is often 
performed by converting the IE into a linear system of equations using the moment 
method (MM) procedure. The MM solution often requires the generation and storage 
of 0(N 2 ) matrix entries and, if a direct solver such as LU decomposition is utilized, 
0(N 3 ) operations are required for its solution, where N is the number of unknowns. 
If an iterative solver such as the biconjugate gradient (BiCG) method is used, the 
solution can be found in 0(r • i • TV 2 ) operations where r is the number of right-hand 
sides and i is the number of iterations per right-hand side. However, for circulant or 
block circulant matrices, the solution may be reached in 0(r * i * N\og 2 N) operations. 
This requires the use of FFTs for computing the matrix-vector products in the BiCG 
algorithm and consequently the resulting solution scheme is often referred to as the 
BiCG-FFT method. In this appendix, we will present specific examples of this 
efficient technique for 2-D and 3-D geometries. We will first look at the relatively 
straightforward 2-D problem followed by the necessarily more involved 3-D case. 

C.l 2-D Integral Equations 

Suppose a flat resistive strip having width w centered at the origin of the y— 0 
plane is excited by an E-polarized (TM Z ) plane wave as shown in Figure C.l. The 
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Figure C.l: Flat resistive strip geometry. 
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appropriate Electric Field Integral Equation (EFIE) may be formed by enforcing the 
resistive transition condition on the strip giving 

E' z (x) = R(x)J t (x) + - J ^ J 2 (x')//< 2 > (k 0 \x - x'|) dx (C.l ) 

where E x z {x) is the incident electric field, R(x) is the normalized resistivity as a func- 
tion of lateral position, J z (x ) is the equivalent electric current on the strip induced by 
the incident field and k Q = y- is the free-space wavenumber. The time dependency 
e 3 * 1 is assumed and suppressed. In (C.l), the primed coordinate denotes the source 
point while the unprimed one indicates the test point. Once the current has been 


determined by solving (C.l), the scattered field in the far-zone is given by 


Ei(4>) 


-Hko P-f ) 



- Z °\l^ J ; JJx')e jk ° T cos * 


(C.2) 


where Z 0 is the free-space intrinsic impedance. 

We proceed with the numerical solution of (C.l) by expanding the unknown 
current in terms of subdomain basis functions as 


TV— 1 


Jz{x) = l 2 J l n }Wn(x ) 


(0.3) 


n=0 


where 


w w 

W n {x) = W(x) nAx — — < x < (n + l)Ax — — 


= 0 


else 


(C.4) 


and Ax = jj. The weight functions VE(x) may assume various forms, one of the 
simplest being a pulse (W(x) — 1). Substituting (C.3) into (C.l) and performing 
Galerkin’s testing, we obtain 

/ i£ TV-1 ✓ i £ 

\W n (x)E[(x)dx = £.7[n'] \ \Wi(x)R(x)dx ,+ 

I t£ m x 

f J-. /_! W4x)W„,(x)HV\K\x - x'l) Jx'dx} 

(C.5) 


n =0 
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which is the discrete form of the integral equation. When pulse basis functions are 
used. (C.5) becomes 

f(n+l)Ar- IJr , f , r(n+l)Ax-f 

I E' z {x)dx = V] J[n K <5[n — n ] / R(x)dx + 

^ JnAx - ~ 

k /<-»>**-¥ 

4 JnAx— 1 ^ J n Ax- y J 

(C.6) 

where the Kronecker delta function 


8[n — n\ = 1 n = n 

f 

0 n 7 C n 


(C.7) 


has been utilized to indicate the self-cell. Assuming that the resistivity is constant 
within a segment (e.g. R[n] = R(x) for nAx < x < (n + l)Az) and making the 
change of variables £ = a: + = x + y, we have 


ytn+ija ^ ^ ^ J[ n '] jtf[n]Ax 6 [n - n'j 

A: r(n+l)Ax Wn' + lJAr , , ^ 

T . 2 ) (*.l« - « IV? 

4 JnAr iii Ai > 


+ 


n =0 


(C. 8 ) 


We now observe that the double integral is in convolutional form and since the 
segments are of uniform length (Aa;), we may introduce the discrete function 


, b f(n+l)Ax An +l)Ax 

s[n - " ] = T / /. H ° ’(‘•I* - 1 !)« 

4 JnAx Jn £ 


’Ax 


(C.9) 


and rewrite (C. 8 ) as 


(n+l)Ax 


r(n+. 

JnAx 


N-l 


AT-1 


£*(£)</£ = Ax Y J[n]R[n]8[n - n] + Y J[n]g[n ~ n ] 

n =0 n— 0 


(C. 10 ) 
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The first sum in (C.10) is recognized as the product of a diagonal matrix ( A.r R{n' Shi- 
rt ]) and the unknown data vector (</[;/]. ;?’ = 0. 1.2 A' — 1) whereas the second 

sum is a truncated, discrete, linear convolution which may be efficiently calculated 
using FFTs. 

The discrete form of the IE (C.10) may be written as a matrix equation 

= {/} (C.ll) 

where the excitation vector entries are given by 

f(n+ l)Ar . / w \ 

/w - L < ci2 > 

and the impedance matrix is given by 

Z[n, n ] = Aa:/?[n]6[n — n ] + g[n — n '] (0.13) 

Alternatively, in preparing to take advantage of the convolution property, we may 
rewrite (C.10) as 

[Aa:#[n]£[n - n']] {J} + [#[n - n']] { J} = {/} (C.14) 

Obviously, the first matrix- vector product can be trivially computed in N operations 

y[n\ = J p [n]AxR[n] n = 0, 1, 2, ..., N — 1 (C.15) 

where the superscript p denotes the iteration. The second matrix-vector product 
involves a fully populated matrix and would thus normally require G(N 2 ) operations 
per iteration for its execution. However, since g[n — n’] is a discrete convolution 
operator, the product may be computed by invoking the discrete convolution theorem 


{ J*} [ 9 [n - n']] = 


(C. 16) 
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where • denotes a Hadamard product and the discrete Fourier transform (DFT) pair 
is given by 

M— 1 

^d{/M} = 

m=0 
i M-l 

{*■[«]} = (C.i7) 

iVI 9=0 

In practice, a radix-2 FFT algorithm is used which requires only C>(A r log 2 A') oper- 
ations as compared to the DFT given in (C.17) which requires 0(N 2 ) operations. 
The pulse function H[N] indicates that only N of a possible M values of the dis- 
crete convolution are retained. Also, the DFT and its inverse operate on periodic 
sequences of period M (all sequences with a tilde are periodic in this appendix). For 
example, the unknown iterate ( J p ) is given by 

j p [n) = J[n] n = 0, 1,2, ..., N — 1 

0 n = N y N + 1, N + 2, ..., M — 1 (C.18) 

It is instructive at this point to examine a couple of fundamental properties of 
DFTs and discrete convolutions. DFTs always assume periodic sequences as shown 
for example in Figure C.2. In this, a single period of the data sequence to be trans- 
formed is indicated by the standard sequence window. This particular data segment 
is the one usually provided to the DFT. However, since the sequence is periodic, any 
full period sequence window is permitted as shown in Figure C.2. Accordingly, in 
the BiCG-FFT method, the user may select which sequence window to use provided 
an appropriate mapping functions between the data sequence and the unknowns are 
specified. This flexibility may ease the programming task as is discussed below. Com- 
putation of a discrete linear convolution via a discrete circular convolution requires 
that M > 2N — 1 values of the matrix sequence be provided. In particular, these 
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standard sequence window 


alternative sequence window 


Figure C.2: Periodic sequence with even symmetry. 
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values must span all possible combinations of source and test unknown interactions. 
The utility of the BiCG-FFT method rests upon the fact that for certain problems, 
a particular redundancy exists within the matrix. For example, suppose the strip 
shown in Figure C.l were divided into N segments numbered consecutively from left 
to right. The resulting matrix will be Toeplitz and there will be N distinct matrix 
entries corresponding to the first row of the matrix. The first row can be loaded into 
the left half of the standard sequence window shown in Figure C.2. The N th row 
of the matrix forms the right half of that sequence and since the matrix is Toeplitz. 
the sequence is even symmetric. Thus, the left half of the sequence corresponds to 
interactions between a test cell and source cells which are physically to the right of 
the test cell. Likewise, the left half of the sequence corresponds to interactions when 
the source cells are to the left of the test cell. The first entry in the sequence is 

the self-cell contribution. Once again, the data sequence contains the full span of 

possible physical alignment between the test and source cells. This example is a case 
where the even symmetry replication rule may be applied 

g[n] = g[n\ n = 0, 1, 2, ..., N - 1 

= 0 n = TV, N + 1, M + 2, ..., — — 1 

= g[M — n] n=y,y + l,|- + 2,....M-l (C.19) 

and we have assumed that g[n — n] = g[n — n]. As shown in (C.19), if M > 2N — lit 
is necessary to pad g[n] with zeros in the middle of the sequence. Use of a replication 
rule reduces the computational burden to the method since only one row of the matrix 
need be computed. For situations where first and last rows of the matrix are not 
mirror images of each other (e.g. it matters whether the source cell is to the right or 
left of the test cell), an asymmetric sequence is requires as shown in Figure C.3. In 
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standard sequence window 

alternative sequence window 


Figure C.3: Asymmetric periodic sequence. 
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this case, all 2.V — 1 interactions need be explicitly computed. 

The 2- D examples given in this appendix all possess even symmetry. However, 
the cross-terms presented in the next section for 3-D problems are not even symmet- 
ric. As it turns out however, these terms do have some redundancy which may be 
exploited by using a non-standard sequence window in order to reduce the matrix 
fill time. The details of such an operation are beyond the scope of this appendix. 

Combining (C.15) and (C.16), the matrix vector product in the iteration cycle of 
the BiCG-FFT method is efficiently computed as 

[Z]{J p } = Ax{J p [n]/?[n]} + //[AT]^ 1 {^ 0 {>}«^{ff}} (C.20) 


This computation requires 0(N + (M)log 2 (M)) operations, provided radix-2 FFTs 
are used in (C.20), and should be compared to the standard matrix vector product 

N - 1 

[. Z]{J P } = £ J p [n']Z[n,n] n = 0, 1,2, ..., N - 1 (C.21) 

n =0 

Figure C.4 illustrates the comparison between (C.21) and the discrete convolution 
(C.20) computed using radix-2 FFTs. Clearly, for N > 30, (C.16) is more efficient 
than (C.21). 

Another 2-D geometry which yields systems which may be converted into circu- 
lant matrices is the circular strip such as the one shown in Figure C.5. The EFIE 
for a resistive circular strip is given by 


E' z (<t>) = R(<f>)J g (<f>) + ^ J ^ J z (<f>')H^ 2) ^2 k 0 a sin ) d 4> d <t> 


(C.22) 


where a is the radius of the arc and a is the subtended angle. Note that if the arc 
is closed, a = 360°. A discrete form of (C.22) may be obtained by using pulse basis 
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Figure C.4: Comparison of operation count for a full matrix- vector product verses a 
FFT-based matrix-vector product for a 2-D formulation. 
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Figure C.5: Circular arc geometry. 
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and uniform zoning (Ao = y). have 


r(r)+l)A<j Q 

/ E\{0--)do 

Jn&<p ~ 


A* - 1 


Y J l u ' 

n'=0 

k 0 a r(*+i)*o n 
4 JnAtp dn* 


/?[n]Ac>^[n — n 

, I 

n =0 

k n a r( n + *)A<> 


//i 2) 2/c.a 


sin 


0-0 


A0 


do do 
(C.23) 


which is similar to (C.8). The entries of the Green’s function matrix are given by 

f , t* > 

0—0 


g[n - n \ 


^ r' B ‘ r+' W H?< 2k»a 

4 Jn&4> Jn &<t> V 


sin 


dxp do 


(0.24) 


Indeed, (C.23) may be solved in exactly the same manner as (C.8) using FFTs. The 
only difference between the flat strip and the circular arc is a possible additional 
symmetry present in the Green’s function sequence. If the arc is closed (a = 360°), 
we find that 

g[N — n] = g[n) n = 1,2,3, ..., y - 1 (C.25) 

which indicates that only y + 1 entries of the sequence need be computed rather 
than N. Therefore, the solution of (C.23) for a — 360° requires roughly half the 
number of matrix entry evaluations as compared to an equal length flat strip. For 
a < 360°, a similar symmetry exists to a lesser extent as long as a > 180°. In this 
case, although the Green’s function sequence (C.24) is periodic, it is incorrect to 
assume that (C.23) now involves a circular convolution. The convolution is still a 
truncated, linear, discrete convolution and in practice the additional symmetry of 
(C.25) is not exploited unless matrix build time is excessive. 
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C.2 3-D Integral Equations 

As was the case for 2-D geometries, there are certain 3-D geometries which admit 
efficient solution methodologies by making use of the FFT-based implementation of 
matrix-vector products. Three of these geometries are the flat plate, an impedance 
insert in a ground plane and an impedance insert in an infinite metallic circular 
cylinder. These geometries are shown in Figure C.6. We shall now proceed with 
the formulation for the later two problems (e.g. the planar and cylindrical inserts) 
using a general coordinate system. In the following, the general coordinates (u,v) 
may be considered as (u = x,v = y) for the planar case and (u = cn£, v = z ) for the 
cylindrical case. 

An integral equation may be obtained by enforcing the standard impedance 
boundary condition (SIBC) 

n x h x E{u,v ) = —r)Z 0 n x H(u,v) (C.26) 


where n is an outward directed unit normal and the total magnetic field is given by 


H(u,v) = H 90 (u,v) + 

jk 0 Y 0 J Gi{u — u, v — v ') ■ [n x E(u , v ) 


dudv (C.27) 


—4 

where S denotes the surface of the insert. In (C.27), the surface field H 9 °(u,v) is 
comprised of the incident field H'(u,v ) and the field reflected by the cylinder when 
no insert is present (H T (u,v)). For the planar geometry, the second kind Green’s 
function is given by 


G 2 {u — u,v — v) = 2G 0 (u — u , v — v ) = 2 


7+ VV 


0 —jk o y/ (u-u') 2 +(v-v ') 2 


k 2 -I 47 T^(u — U ) 2 + (v — v ) 2 

(C.28) 
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Figure C. 6: 3-D geometries: (a) flat impedance insert in a groundplane (or 

impedance plate in free-space) (b) impedance insert in a circular cylinder 
(2 a = A <t>). 
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where G 0 is the free-space dyadic Green's function and / = xx + yy + zz is the idem 
factor. For a cylinder, the appropriate Green's function is available as an eigenvalue 
series [29] for small radius cylinders or as a creeping wave series [36] for large radius 
cylinders. Irrespective of the form of G 2 , upon making use of (C.27) into (C.26). we 
have 


^ * . riao, ^ nx fix E(u,v) 
-Z 0 n x [U'V) = h 


jk 0 fl X / G2 (u — u , V — V ) * 77 X E(u , V ) 

Js 


du dv 


(C.29) 


We note that for planar geometries, an appropriate right-handed system is (u,t), n) 
whereas for cylindrical problems, the right-handed system is given by Upon 

decomposing the fields (E,H 9 °) and the dyadic Green’s function (G 2 ) in terms of 
tangential components (it, v), (C.29) yields the following coupled set of integral equa- 
tions 


u:Z 0 H*°{u,v) = 


E u (u,v) 


+ 


t 

jk 0 J s ^E u (u ,v')G uu (u -v! ,v - v) - 


v : -Z 0 H 9 u °(u,v) = 


E v (u ,v')G uv (u — u , v — v ) 
E u (u,v) 


du dv 


+ 


jk 0 v')G vu (u - u, v - v') - 

E v (u ,v')G vv (u — u ,v — v) 


du dv (C.30) 


To convert (C.30) into a discrete set of equations, we expand the unknown electric 
field components in terms of subdomain basis functions 

N —3 N — 2 

E u (u,v ) = jr jh E u [t,s]W^ 3 (u,v) 

t = 0 5=0 
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E v (u,v) = E v [t.s]\\?{u.v) ( c:m ) 

f = 0 5 = 0 

where N u denotes the number of nodes in the u-direction and A’,, is the number of 
nodes in the n-direction. As shown in Figure C.7, t denotes the row number of t ho 
edge discretization and s is the column number. For this example, the unknowns are 
associated with the free-edges and the associated basis functions may be represented 
as 


W^ s (u,v) = j if local edge = 1 

(v t - v b ) 

= 7 —^ \ if local edge — 2 

{v t - ut) 

0 else 


«'>,») = 


(u r — u) 
(u r - Ul ) 
(tt - Ul) 
(Mr - Ul) 


i f local edge = 3 
if local edge = 4 


0 else 


(C.32) 


where the local edge numbering is illustrated in Figure C. 8 . Note that a free edge 
is one that is not tangential to any metallic walls and that use of (C.32) requires 
a finite element type assembly. For the discretization shown in Figure C. 8 , there 
are ( N u — l)(N v — 1 ) elements and a total of 2N U N V — 3(N U + N v ) unknowns. Of 
these, ( N u — l)(Af„ — 2 ) unknowns are associated with the u-directed field while 
(N u — 2)(N V — 1 ) unknowns represent the u-directed field. 

Substituting the field expansions (C.31) into the coupled integral equations (C.30) 
and employing the Galerkin’s testing procedure, we have 

N v -3N u -2 

F u [t,s] = “ 5 ] + 

t'= 0 s'=0 

N v — 2 Nu - 3 

Y, E E v [t\s}g uv [t-t\s-s] t = 0,l,2,...,fV„-3 

t f =0 s' =0 
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* = 0.1,2 Y u - 2 

Sr —3 Su —2 

F„[f.s] = ^ ^ E u [i .s]g vu [t - t - s] + 

t'= 0 a'=0 
Sr— 2 A’ u -3 

53 E -t',s - s'] t = 0,1.2 v,, - 2 

('=0 a '=0 

5 = 0, 1,2, ..., A„ - 3 


(C.33) 


where 


Quu t ,s s ] 

SW t , s 5 ] 

t r S S ] 
Qvu t ) S S ] 

F U [M] 

F„[M] 


f W u (u ,v')W u {u,v) 
JS' 


, dudv 


u = u 
1 

v — t> 


6 [e — e 

rj[e) 

—jk 0 f [ W u (u ,v')W u (u,v)G uv (u — u ,v — v )du dv dudv 

Js e Js , 

6[e - e] 


r)[e] 


( W v {u ,v')W v {u,v) 
Js c 


, dudv 


U — u 
/ 

V — V 


+jk Q f / WJu ,v')W v (u,v)G vu (u - u,v - v')du ' dv 'dudv 
Js e JSj 

jk Q f [ W v {u ,v')W u {u,v)G uu (u - u',v - v')du 'dv 'dudv 

Js e Js el 

jk 0 f f W u (u ,v')W v (u,v)G vv (u - u',v - v')du 'dv 'dudv 
Js c Js t , 

Z 0 [ W u {u,v)H°°(u,v)dudv 
J S a 

-Z 0 [ W v {u,v)H 9 u °{u,v)dudv (C.34) 

J Sc 


and e refers to the test element while e denotes the source element. We note that 
each of the double sums in (C.33) is a truncated, discrete, linear convolution and 
hence amenable to the BiCG-FFT method. 

To proceed, we define the 2-D DFT pair (analogous to (C.17)) 


^2D{/M} 
?2D {F[<?,P]} 


-(sp+tq) 


Mj-1 M 2 -l 

E E f[t,s]e->^ 

t - 0 s - 0 

1 Ml— 1 M 2 — 1 . 2 

— — Y- y F[q,p]e 3 -^ {3p+tq) 

Ml M 2 £5 l9,PJ 


(C.35) 
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Using (C.35). the convolutions in (C.33) can then be rewritten as 

Mi — 1 M 2 -l 

E E C[lV] 9 [f {E{t.s}} 'F,n {»[(..<]}} 

t’= 0 s=0 

(C .36) 

The order of the relevant DFTs must be M\ > 2(number of rows) - 1 and M 2 > 
2(number of columns) - 1 where the number of rows and columns of the discretization 
may vary with each convolution (see Figure C.7). For example, the first convolution 
in (C.33) is associated with {/-testing and i/-source edges and hence the number of 
rows and columns is (N v — 2) and ( N u — 1), respectively. The field sequences are 
loaded into a M\ x M 2 array in row/column order of the field discretization and 
the remaining entries form a zero pad. The Green’s function sequence must be 
loaded into a similar array (in the same manner) and periodic replication must be 
performed to provide the necessary “negative lags”. If the sequence has the property, 
g[t — t\ s — s] = g[t‘ — t,s — 5 ], then this replication process takes the form 

g[t, s ] = g[t , s] 0 <t <— - 1 0 < s < — -1 

i 

= g[M\ + 2 -t,s] y^K^-lO^Ky-l 

= g[t, M2 + 2-5] 0 < f < — --l-y<s<M 2 -l 

Mi jt/ 1) 

= + 2 - M 2 + 2 - s] — - < t < M 1 - 1 -y < s < M 2 - 1 

(C.37) 

For the presented example, g uu [t, s] and g vv \t,s] possess this property while the cross- 
terms, 5r U t)[U' s ] and fi'vuiU 5 ]? may not. In the case of a planar insert, the cross-terms 
do possess even symmetry due to the use of a symmetric closed form dyadic Green’s 
function. On the other hand, the cylindrical insert does not have such symmetry since 
the equidistant interactions may differ by a sign on account of the form of the dyadic 
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Green's function. For example, relative to the first quadrant of the array, interactions 
in the second and fourth quadrants have the opposite sign while the third quadrant 
have the like sign. These signs may be accounted for by using a complex replication 
rule and an non-standard sequence window. Otherwise, all possible lags must be 
computed requiring longer matrix build time since more than the first u-directed 
and e-directed edges need be used as sources. 

Once the periodic arrays are loaded, the required matrix-vector product for the 
uu-interactions may be performed in 0 ((MilogMi)(A/ 2 logA/ 2 )) operations rather 
than the 0(((N U — 2)(N V — 3)) 2 ) operations required for a standard matrix-vector 
product. The comparison is shown in Figure C.9 with Mi = 2 (N v — 3), M 2 = 
2(N U — 2) and N v ~ N u = N. Clearly, when the number of nodes per side exceeds 
10-15, the FFT-based matrix-vector product is more efficient than a conventional 
matrix-vector product. In practice, the FFT-based product is more efficient than a 
standard product in terms of wall clock time for N < 10 since in order to exploit the 
memory savings afforded by uniform zoning of a convolutional kernel without using 
FFTs, additional overhead is incurred to match the appropriate matrix entry with 
the correct vector entries. Similar results are obtained for the other convolutions in 
(C.33). 
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Number of Rows/Columns in Grid (N) 


Figure C. 9: Comparison of operation count for a full matrix-vector product verses 
a FFT-based matrix-vector product for a 3-D formulation. N is the 
number of nodes in each direction of the grid, M\ = 2(N-3) and M 2 = 
2(N-2). 
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